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I. INTRODUCTION 


Radome design is a composite discipline evolved from 
aeronautical, electrical, mechanical and materials engineering. The 
word radome was coined during World War Il as a composite of two 
words, radar and dome [Ref. l:p. ix]. The term radome originally 
described the family of dome-shaped structures designed and built 
to house and protect the radar antennas on aircraft, but now the 
term is applied to any structure that protects a radio frequency 
antenna system. The structure is designed to provide maximum 
protection. to the antenna with minimum interference to the 
antenna's electrical performance. 

Radomes on missiles are presently designed to meet many 
diverse specifications. These specification requirements can be 
divided into two categories, electrical and aero-mechanical. Electrical 
design problems include those that have to do with the radome's 
ability to act as a window, ideally transparent, to the electromagnetic 
energy that must pass through the radome to be used by the 
missile's antenna system. This energy is of a specific frequency 
range and intensity that is within the operating parameters of the 
antenna system. The magnitude and form of the electromagnetic 
energy will be dependent on range from the source, polarization, and 
incident angle. Aero-mechanical design problems address not only 


the design of the structure of the radome, but also the ability of the 


radome to withstand environmental conditions and battle damage. 
Environmental conditions vary with scenario, and include: normal 
pressure force and viscosity, wind loads, inertial loads, temperature 
variations, moisture resistance, resistance to rain erosion, and other 
diverse factors. Design of a radome for use in battle may require 
that the aero-mechanical design consider resistance to nuclear 
radiation and laser weapons. 

As in any engineering design, there are tradeoffs between the 
various requirements, such that each interrelated component design 
factor must be weighted and considered in the overall system design. 
Numerical models of environmental conditions assist in the analysis 
of the tradeoffs for the aero-mechanical design. Numerical models to 
aid in the electrical design analysis of a radome to protect missile 
antenna systems which operate at high frequencies are available, 
such as the ray-tracing algorithms described by Hirsch and Grove in 
Ref. 2. When a missile antenna system requires operation in the 
resonance region of frequencies for the radome parameters, a 
numerical model is required that can solve the radome interaction 
problem. 

The subject of this thesis is the development of a generalized, 
interactive set of programs to be used in the analysis of the 
interaction between an actual or proposed radome structure and 
incident electromagnetic fields at the lower frequencies of interest, 
where strong resonance characteristics exist. In order to analyze a 
specified radome, a numerical model is constructed using the 


interactive program, RADSHP. This radome generation program 
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provides geometric shaping functions to allow the user to construct a 
multilayer radome model of specified materials. Chapter II provides 
a discussion of radome design methodology, the aerodynamic shaping 
functions available in this package, radome wall constructions, and 
material specification. 

This thesis work is an application of the programs developed 
for use on a personal computer by Morgan [Ref. 3] and Connolly 
[Ref. 4] to calculate electromagnetic scattering by inhomogeneous 
dielectric and/or magnetic bodies of revolution. As discussed in 
Chapter Il, radomes are manufactured using dielectric materials in 
order to provide a radio frequency window to the antenna housed in 
the radome's interior. An axiSymmetric missile radome structure 
meets the physical structure requirement for rotational symmetry. 
The numerical analysis of the interaction of the specified radome 
structure and the incident electromagnetic energy 1s performed by 
the code developed by Morgan, EMCAD. Modifications to this code 
resulted in the program EMRAD, which provides data file generation 
to store the modal expansion coefficients which are used by the 
program CORFLD to determine the electromagnetic fields within the 
radome interior. In Chapter III, the theory required for the 
formulation of the computed interior fields of the radome by the 
program CORFLD is discussed. 

The programs developed by Morgan use an optimized 
variational finite-element algorithm in conjunction with a 
tri-regional unimoment method to provide scattering solutions for 


each of multiple incident fields impinging on the specified structure. 
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The programs represent a significant improvement in computational 
speed and versatility over computer codes which find the solution to 
the electromagnetic resonance behavior by use of surface current 
integral equation formulations. An overview of the finite element 
scattering algorithm, along with the concepts embodied in the 
coupled azimuthal potentials and the unimoment method are 
included in Appendix A. Further information is available in 
references 3, 4, 5, 6, and 7. 

Chapter IV discusses the flow of data between the programs 
which have been developed. The general use of the RADSHP 
program to construct the radome data files is discussed, as well as 
the data files generated and their uses. The flow of information into 
and out of the radome analysis program EMRAD is analyzed. The 
execution and operation of CORFLD is discussed. 

In Chapter V, the validation of the interior field computation is 
discussed, wherein the numerically generated fields are compared to 
the incident fields for a radome constructed of air. The results 
produced for radomes constructed of materials presently used in the 
radome industry are presented and discussed. 

The conclusions in Chapter VI will summarize the work 
previously presented, as well as provide recommendations for future 
modifications to the interactive package of programs being used. 

Additional Appendices are included to provide copies of the 
source code for the computer programs, RADSHP and CORFLD, which 
were written as part of this thesis effort. Source listings for 
additional programs, such as EMCADIN, EMESH, and EMRAD are 
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available upon request from Prof. Michael A. Morgan. His address is 


given in Appendix D, Software Sources. 


II. RADOME DESIGN 


A. RADOME DESIGN METHODOLOGY 

The configuration of a radome for high-speed missiles must 
combine high aerodynamic performance with a highly accurate 
transmission of electromagnetic energy through the radome to the 
missile seeker. As described in Ref. 2, the methodology for 
establishing radome performance requirements is an iterative 
process. Figure 2.1 shows the process as a set of design procedure 
loops that begin with a statement of desired performance, and 
include analytical studies, hardware development, and hardware 
evaluation. These result in a set of optimized radome requirements, 
which is the blueprint from which a radome can be constructed. 
Computer software tools help to maximize the benefit of the work 
done in analytical studies, thereby minimizing the number of 
iterations required for the hardware development and evaluation. 
The programs described in this thesis are designed to analyze the 
interaction between the aerodynamic considerations in radome 
design and the electrical performance of the radome. 

The aerodynamic considerations in the design of a radome 
Structure are three-fold: shape-generated drag, thermal stress from 


aerodynamic heating, and structural performance. These three 
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Methodology for establishing radome requirements 
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parameters are not independent of each other, and the electrical 
performance of the radome is a strong function of these aerodynamic 
factors. The outer shape of the radome has a large impact on how 
the missile will fly through the air. It is this outer conformation that 
will, along with speed of travel and angle of attack, determine the 
forces that act on the radome structure during flight. The forces 
generated during flight determine the design requirements on the 
radome structure for material strength and thermal stress. 
Fortunately, the shapes that result in the best drag performance are 
those that can also be used to develop radomes that match the 
structural requirements. Unfortunately, the shapes that best meet 
the drag requirements are not necessarily the best for maximizing 
the electrical performance. The programs developed in this thesis 
include the computer aided design (CAD) functions to analyze several 
Shaping functions that aerodynamicists have developed to minimize 
the aerodynamic drag forces with respect to radome structure 


parameters. 


B. RADOME SHAPING FUNCTIONS 

An ogive shaping function is one of the most common used in 
radome designs. The ogive shape is popular because it is easily 
fabricated and has relatively acceptable aerodynamic performance at 
high speeds. It is defined [Ref. 1:p. 48] as a segment of the arc of a 


circle, whose radius, F,, is greater than the radius of the radome 
base. A pointed ogive is generated when the circle's radius, F,, is 
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large enough so that the arc reaches the radome's axis of symmetry. 


The ogive shaping function is described by 


x= Fi -(H-)-F} -F . (2.1) 
where F, and F} determine the position of the circular shaping arc, as 
illustrated in Fig. 2.2. The ogive shaping function degenerates to a 
hemisphere when the shaping factors, F, and F,, are defined as 
F, =F, and F, = 0. 

A tangent ogive shaping function generates a radome shape 
such that the surface of the radome is tangent at the base of the 


radome, providing a smooth transition to the missile body [Ref. 8]. A 


tangent ogive shaping function is determined such that 


x =4 (Fez = Fs | (2.2) 


where 
R H? 
Fy = — 2 — 25 
4 > R >) 
and 
Fs, = Fy - - ! (2.4) 


A shaping function that provides maximum  volume-to-drag 


ratio at high flight speeds is the von Karman function [Ref. l:p. 48] , 


Nye 





C Im 


l lewis 
RJ E c 7 sinQE) | ; PAS) 


9 


ET] 





Figure 2.2 Radome shape - nomenclature 


where 


BE = cos) = , E2 | (DEO 


Radomes described by this von Karman shaping function have a 
limited application due to fabrication difficulties related to the 
complex shape. 

In order to minimize drag for a given height-to-radius ratio, a 
Newtonian contour is used. This contour is approximated by the 


power series radome shaping function [Ref. 1:p. 49], 


F6 
XE R(E? | , (2.7) 


The exponent, Fg, provides for shaping. When Fẹ = 1.0, the radome 
shape generated is conical.  Aerodynamicists normally use F¿ = 0.75 
to approximate the Newtonian contour. 


A shaping function analyzed by Hirsch and Grove [Ref. 2] uses a 


parabolic basis function for the radome shape, 


x? = F,-z. (2.8) 
The rounded effect developed at the tip of the radome by the 
parabolic shape function can be eliminated by the creation of a 
mathematical "deadzone" for points within a region about the 
radome's apex. The parabolic shaping coefficient, F,, determines the 
extent of the "dead zone" and thereby controls the slope of the 
parabolic surface. The radome shape is determined in terms of the 


shaping coefficient, then scaled to the radome size parameters. The 
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unscaled shaping function, 


X = JF; LN (2.9) 


is used to generate a parabolic radome shape with a height of F,. In 


order to eliminate the parabolic rounding of the apex, the “deadzone" 


is established for any points for the unscaled z value within 1.0 of F;. 
Within this deadzone, z > ( F,- 1 ), the radome shape is determined 
at a line of constant slope, converging at a height of z = F,+ 1.0. This 


results in the marrying of a pointed apex onto the parabolic radome 


Shape, with the magnitude of the shaping coefficient, F,, determining 


the relative magnitude of the apex to the total radome height. 

This project contains coding in RADSHP to generate radomes of 
multiple layers for the shaping functions described: the general 
ogive, the tangent ogive, the von Karman shape, the power series 
radome, and the parabolic radome with pointed apex. RADSHP can 
generate a complex radome construction model using different 
shaping functions for each layer. RADSHP generates a material data 
file appropriate for input to the numerical analysis program, EMRAD, 
and a structure data file that can be viewed and modified using 
CURVE DIGITIZER (Appendix 2) , then translated by EMCADIN into an 
input structure data file for EMRAD. 


C. RADOME WALL CONSTRUCTION 
If a material existed which was mechanically strong enough to 


sustain the aerodynamic forces that the radome experiences during 
eZ 


flight, and that same material had the electromagnetic properties of a 
vacuum, then the radome could be built to ideally meet all 
requirements. “Such a material does not exist, and, in general, as the 
requirement for material strength increases, the electromagnetic 
properties tend to strongly deviate from the properties of the ideal 
vacuum. The transmissivity of the radome tends to decrease as the 
dielectric constants of the construction materials increase. In order 
to meet strength requirements, while minimizing the perturbation 
caused to the electromagnetic field, layer construction of radomes 
has been developed. 

Lamination construction methods allow for an alternating 
combination of thin layers of high density, high strength materials of 
high dielectric constants, with thicker layers of low density, lower 
strength materials of low dielectric constants. Electromagnetic design 
of laminated radomes is more complex than the design of monolithic 
radomes in that there is a greater number of boundary conditions to 
be considered in establishing the lamination layer thicknesses. The 
EMRAD program developed by Morgan [Ref. 3] currently handles 
layered dielectrics with up to five layers. 

Figure 2.3, adapted from the Radome Engineering Handbook 
[Ref. 1:p. 46], shows several types of laminated wall construction that 
exhibit the alternation of layers of low dielectric constant and high 
dielectric constant materials. The two-ply sandwich is the simplest 
form of lamination construction with a thin skin of strong, dense 


material of high dielectric constant, supported internally by a thicker 
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layer of porous material of low dielectric constant. This design is 
especially attractive at low frequencies as the lateral dimension of 
the thin layer of dense material is only a fraction of the longer 
wavelength. 

The A-sandwich and B-sandwich are two variations of a 
stronger three-layer construction. Although the A-sandwich is 
Superior electrically to the B-sandwich, while the placement of the 
high density material in the center of the layer construction will 
Strengthen the radome wall by making it stiffer, the porous, low 
density material on the exterior of the radome wall will not sustain 
the aerodynamic forces during flight. 

The B-sandwich's inner and outer skin of high density material 
and the core of low dielectric material provide stiffness and strength 
to the radome structure while providing a dense outer skin. This 
construction has a higher strength-to-weight ratio than a monolithic 
construction. The core thickness is determined by the strength 
requirements and provides spacing between the inner and outer 
Skin. Optimal electrical spacing between the skins would by 
frequency dependent at a core thickness of a quarter wavelength 
[Ref 1:p. 45]. 

The C-sandwich and multilayer radome wall constructions 


provide for very high strength levels which may be required at high 
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Figure 2.3 Laminated radome wall construction 
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Operating speeds or for use in inclement weather. The layer material 
may consist of plastic laminates and ceramics. Multilayered 
construction is attractive where extreme structural rigidity or 
broadband capabilities are primary requirements. [Ref. 1:p. 46] 

The radome interacts with the electromagnetic field as a lens, 
modifying the field as it passes through the radome structure. At 
higher frequencies, this interaction can predicted by the use of 
ray-tracing techniques, similar to optical lens analysis [Refs. 2 and 
8]. Studies have been conducted at higher frequencies to design a 
layer of dielectric material within the radome's interior structure to 
act aS a correcting lens [Ref. 8] as shown in Fig. 2.4 . This package of 
programs will allow similar analysis of lens construction at lower 
frequencies, where the radome and lens parameters are within the 


resonance region of the electromagnetic fields. 


D. MATERIALS FOR RADOME CONSTRUCTION 

The choice of material for construction of a radome is another 
engineering design decision. As material science engineers continue 
to develop new materials, particularly low-weight, high-strength 
ceramics, the range of materials available to the radome designer 
increases. The general properties of the three groups of materials 
which are in common use are listed in Ref. 9, and is used here as a 
source for material selection. The groups of materials discussed are 
fiber-reinforced organic resins, refractory oxides, and glass ceramics. 


Physical properties are listed in Table 2.1. 
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Figure 2.4 Interior radome lens of dielectric material 
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Relative 
permittivity 
(E) 


Loss tangent 
( tan 6 ) 


Density 


Ultimate 

Tensile 

Strength 
(UTS) 


Compressive 
Strength 


Maximum 
Temperature 


TABLE 2.1 


RADOME MATERIAL PROPERTIES 


FIBER 
REINFORCED GLASS REFRACTORY 
RESINS CERAMICS OXIDES 
4-6 5:5 8-9 
0.01 0.003 0.002 
N/A 160 1b/ft 220 Ib/ft? 
40 kpsi 7 kpsi 27 kpsi 
20-40 kpsi 35 kpsi 280 kpsi 
260°C 500 ° C 1700 ° C 
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In fiber-reinforced organic resins, a variety of resin types are 
used to impregnate fiberglass. Organic resins, such as phenolic, 
polyester, and epoxy are used with the fiber to construct the radome. 
They are then thermally treated to set the structure. Several 
manufacturing techniques are discussed in Ref. 1. Control of layer 
thickness is difficult to maintain and requires finishing to avoid 
aberration. Radomes constructed of these materials have poor 
resistance to rain erosion, but are relatively cheap and easy to 
manufacture, even in large sizes. Maximum temperatures that the 
material can sustain limit use of these materials to missiles with 
speed of less than Mach 1.5 to Mach 2.0. 

Typical glass ceramics used in radome construction are 
Pyrocream and glass-bonded mica. Compared to organic radome 
materials, ceramic materials have much higher resistance to rain 
erosion at supersonic speeds. Ceramics are brittle materials and the 
controls needed during manufacturing to ensure structural integrity 
of the radomes makes them expensive to build [Ref. l:p. 242]. 

The refractory oxide used in radome construction is usually 


alumina (AL,0O4), but magnesia (MgO) and fused silica (5105) are also 


used. Refractory oxide manufacturing is also expensive due to the 
high temperatures used in forming these materials to the desired 
radome shape. These materials have high dielectric constants, that 
can undergo a change of up to ten percent between room 
temperature and 1700 °C. Some of the methods by which refractory 


oxide radomes are manufactured are slip casting, isostatic pressing, 
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and flame spraying [Ref. l:pp. 298-313]. Oxides have a high 
resistance to rain erosion and wall thickness can be controlled to 
within 0.002 inches, allowing for extreme control of aberration. 
Although the dielectric constant of these materials is high, they 
provide the highest degree of thermal protection and place the least 


restriction on missile speed. 
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III. ELECTROMAGNETIC FIELD THEORY 


A. BACKGROUND 

The analysis of the interaction between a missile's radome and 
the electromagnetic signal that is radiated or received by that 
missile's radar is a problem similar to those undergoing intensive 
investigation in electromagnetic interaction with other physical 
medium. Ongoing study of medium interactions include other 
military weapons and weapons platforms, forested areas, bodies of 
water, rain, and biological tissues [Refs. 3, 5, and 10]. The study of 
electromagnetic wave penetration through missile radomes allows 
for a special approach to the solution of the problem for the common 
case of rotational symmetry in the construction of the radome. In 
the high-frequency regime, the solution may be found as described 
by Hirsch and Grove [Ref. 2], with application of an asymptotic 
approach using the geometrical theory of diffraction (GTD), which 
incorporates ray-tracing. When lower frequencies are considered, 
and the interaction problem is in the resonance region, a complete 
solution of Maxwell's equations is required. The asymptotic GTD 
method yields unsatisfactory results. 

The solution of Maxwell's equations for the radome used within 
this project is divided into three regions as illustrated in Fig. 3.1. The 
three regions are the homogenous spherical inner radome core 
(designated Region I); the transition region, where a semi-annular 
conformal mesh is generated for field analysis (designated Region Il); 
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Figure 3.1 Three spatial radome regions for 
electromagnetic field analysis 
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and an unbounded exterior scattering region (designated Region III). 
The total solutions for the fields in the core, Region I, and the 
scattered fields in the exterior region, Region III, are produced via 
separate spherical harmonic expansions. Using the unimoment 
method (Appendix A), the modes in these expansions are applied as 
boundary conditions on the inner and outer spherical surfaces 
bounding the mesh and the resultant fields for each such mode are 
found using the finite element method. A similar analysis is 
performed for the incident fields on the outer surface. 

Using superposition, the finite element solutions are then 
combined using the same unknown coefficients as in the original 
expressions. The combined finite element solution is equated, in the 
least squares sense, to the original expressions along the inner and 
outer surfaces. Enforcement of this equality allows solution of the 
expansion coefficients for the core and exterior regions. 

The exterior expansion coefficients are used by EMRAD to 
assemble the scattered fields and bistatic radar cross sections (see 
Refs. 3 and 4). The core region spherical harmonic expansion 
coefficients are used to construct the fields that penetrate the 
radome to illuminate the enclosed antenna. The details of how these 
interior region fields are generated, using the core region spherical 
harmonics, is the topic of the next section. An overview of the 
theory used by Morgan in EMRAD is included in Appendix A. 


Further information is available in Refs. 3, 4, 5, 6, and 7. 


23 


B. GENERATION OF THE CORE REGION FIELDS 

In this section we will consider the construction of the core 
region fields using the expansion coefficients provided by EMRAD. 
These fields are assembled and plotted in various ways over 
user-selected planar surfaces by CORFLD. 

Within the spherical core, the source-free phasor Maxwell's 


curl equations, with suppressed eJ®! time variation, are 


—-VxE=jw y H (308 


and 


VxH-joeE , (3:25) 
where u and & are the respective permeability and permittivity. The 
u and € will usually be those of free-space for the case of the core 
region in a radome. 

In a homogenous, source-free spherical region, such as the 
radome core, it can be shown [Ref. 11] that the vector 
electromagnetic field can be represented by the combination of two 
types of spherical harmonic fields, transverse magnetic (TM) and 
transverse electric (TE). Using the spherical coordinates of Fig. 3.2, 


these TM and TE fields have the respective properties that H, = O and 
E. = 0. Furthermore, the TM and TE vector fields can be generated 


r 


from radially directed vector potentials, such that for the TM fields, 


A =A,r (3.3) 


and 
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Figure 3.2 Spherical coordinates with unit vectors 
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F = b pe (3.4) 
where the radial components satisfy 
v2 (At) He (4) = 0 (3.5) 


and 


y? E E E ND (3.6) 
T T 


These equations are in the form of the complex scalar 


Helmholtz equation, 


V2w 4k?w-0. (3.7) 


In spherical coordinates, the Helmholtz equation becomes [Ref. 11], 


l 9 6 |. 9 s er PY O . (3.8) 


==» = — n — 
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To find a solution to this differential equation, the method of 


separation of variables is used. Letting 


Y = Rr) 0) O@) , (3.9) 
the separation of variables procedure yields the system of ordinary 


differential equations: 


d | 5 dR DN a. 
af A T -n (nH) |R =0, (3.1 0) 
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lo df. de m? 
Fine ¿sino nan ale ==, (3-1 1) 


and 


——im*® = 0. (3.1 2) 


The R equation is related to Bessel's equation, and its solutions are in 


the form of spherical Bessel functions. Spherical Bessel functions, b, , 
are related to ordinary Bessel functions, B,, as indicated. 


RG) - b,(kn » |—— B 


JF p : (3.13) 


The © equation is related to Legendre's equation and its solutions are 


the associated Legendre functions, 


G(0) = L (cose) > P"(cost) , Q” (cos0) , (3.14) 
where P,m(cos 0) are the associated Legendre functions of the first 
kind and Q,™(cos 80) are the associated Legendre functions of the 


second kind. The associated Legendre functions of the second kind 
tend toward infinity at points along the z-axis and must be excluded 


for solutions which include this region. 


2 


The equation has harmonic sinusoidal type solutions 


D(H) = h(md¢)—> cos(mo) sin(mo), e?/"^ 


c; cos(m $)* c5sin(m Q) 


= Cen A (3.15) 
Together, these separated solutions form the product solutions 


to the Helmholtz equation as 


Tien = Dn ND P^" (cos0) h(md) . (3.1 6) 
Equations of this type form the elementary spherical harmonic wave 
functions. The general solution to the Helmholtz equation is formed 
by summation of linear combinations of these elementary functions, 
one for the TM fields and one for the TE fields. Each of these 
solutions will be of the form 


Y(r,0,9) = Y Y b, (kr) P”(cos8) h(mg) . (3.17) 


m n 


The summation over all possible m and n is limited by the 
nature of the Legendre polynomial of the first kind in that P," =0 
for n < m. The summation is also limited by the restriction that n 
will not equal zero when m is zero. For both n and m equal to 
zero, then the field evaluated is zero, therefore this solution is also 
excluded. Then m is zero and the summation over n will start 


with n equal to one. This leads to the scalar field expansion 
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Y (1,0,0)= Y S Cmn Cosmo +d y. sinmp) P% (cos) b, (kT). (3.1 8) 
m=0 n=m 
nO 


As the Helmholtz equation under consideration is given as 


V^w +k? Y =0, (3.1 9) 
then for the TM case 
y= I l (3.20) 
and for the TE case 
T - i (3.2 1) 


The spherical harmonic solutions of the Helmholtz equations for 
these sets of modes are then of the form 
Ar 
= Tan = [(kr) b,(kD] P,'(cosé) h(mo) . (5092) 
F, 
In this renormalized form of the solution to the scalar Helmholtz 
equation, the Riccati (or Schelkunoff) spherical Bessel functions are 


present. The Riccati spherical Bessel functions relate to the spherical 


Bessel function and normal Bessel functions as follows: 


1 (kn = [(kr)b, (kr)] = fake à, je (3.2.3) 
2 


Therefore, in the spherical core region, the solutions for A, and F, are 
of the form 
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Ar oo. E 

=) »(Cmncosmó +0, sinmg) PP(cos0) J, (kr) . (3.24) 
F m=0 n=m 
: nz 


When expressed in terms of exponential Fourier series in 0, the 


solutions are of the form 


A, oo oo 

=> D ann PACo EA (3.2.5) 
E In--ee nm 
: nz0 


The coefficients within the two summations are determined by the 
solution of the scattering problem provided by EMRAD. 

The spherical coordinate electric and magnetic field 
components can now be found for the TE and TM cases. For the TM 


case, where 


Hrm = Vx(A,r) (3.26) 


and 


pu xo; VXVX(A r), (3.27) 


the field components are 











1 9? 2 
A 
a 31229 
Pom joe r oro0 
PASA 
E eee r 3-30 
IM joe rsin8 9rd6 f ) 
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"TM 

__1_2A, 
Hory - rsin8 00 
a 

E 1099. 


For the TE case, where 
HTE = — Vx( F r) 


and 


] 
= — VxVx(Fr) , 
Erg jue x x( R ) 


the field components are 


| OF 





"OE: 
TE jog r drae 
MORE ee ee, 

TE jmp rsin® drop ` 








The total fields are the sum of the TE and TM parts 


$1 


(23 1) 


(3.3 2) 


(S333) 


(3.34) 


(3.3 5) 


(3.3 6) 


(65377 ) 


(3.3 8) 


ro) 


(3.4 0) 


(3.4 1) 























B= E A (3.42) 

H, = roi E en, (3.45) 
A o 
me d Ue wl INNO (3.47) 


r jou rsin8 oro — 


The solution for the coefficients in the A, and F, expansions is 


performed by the finite element and unimoment solution method 
explained in Appendix A. The finite element method is used by 
EMRAD and applied to Region II between the inner core and outer 
exterior region. In this finite element region there is the 
inhomogeneous axially symmetric structure of the radome. 

EMRAD generates solutions for incident plane waves with 
Poynting vectors at incident angles specified by the program user, 
such that the Poynting vector's direction is specified by the unit 
propagation vector p =sinax +cos az. As any plane wave can be 
decomposed into the sum of two orthogonal, linearly polarized plane 
waves, EMRAD evaluates two plane waves for each incident angle, a 
TM plane wave and a TE plane wave, as shown in Fig. 3.3. The TM 
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Figure 3.3 Incident Poynting vector with TE and TM plane waves 
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and TE plane waves are orthogonal to the Poynting vector and to 
each other. The TE plane wave 1s always in the y direction. The TM 
wave is in the x-z plane and directed as determined by the 
expression cos ax -sina z. For incident plane waves specified in 
this manner, the solution for the component values of the E fields in 
the core of the radome can be determined from taylored equations 
for the TM and the TE incident plane waves. For the TM incident 


wave, the core E-field component values are determined as: 





_ —jE,¿sin0 Pr" (cos0) 
E r8 D as $ acoso È os n (n +)J (kr Dc | (3.4 8) 
nx*0 
Er, - ed Fafe Y e cosmo Y. [ad jo ij dcr Eo 49) 
nz0 


and 


Eyr, = Er Rum S cqsinimy Sa, J kK ab, i (k p) —ų——— PES | (3.50) 
nz 
For the TE incident wave, the core E-field component values are 


determined as: 


M 
ri yy A sino X [s nee ee eee rat Aqu (3.51) 
nx*0 
o d p" 
E,(r,6)- 2L Efe Y 009 3 ad LEE n tb, Jy (k )—72- lo 52) 


nx*0 


and 
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E,(r, TET $ m bosco a ¡Ia ab melas | 
i 


(32559 


With the a, and b, coefficients supplied by EMRAD for each pair of 


plane waves, TM and TE, at each user-specified angle, CORFLD 
generates the interior phasor field components using the equations 
above. The total magnitude of the complex vector field at any point 
within the radome's core is expressed as the square root of the sum 
of the squares of the complex absolute values of the field 
components. This total magnitude is computed and displayed for 
each incident wave by CORFLD, using the a, and b, coefficients 
supplied by EMRAD. 

In order to determine the perturbation of the electromagnetic 
field due to the presence of the radome, CORFLD generates the values 
of the incident plane wave, then compares this to the core region 
field components calculated by EMRAD and CORFLD. The formulas 


for the undisturbed incident TM and TE planes waves are 
EM-E eJKoOPT (xcosa—zsina) (3.54) 
and 
EE-E, eSP y (3.55) 
The magnitude factor, E,= 1.0, is set by EMRAD for both TM and TE 


waves. The unit propagation vector, p, includes the direction of 


propagation for the incident waves and is designated as 


= (x sSina+z cosa) . (3.56) 
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This vector 1s dotted with the position vector, r, for each point on the 


plane of interest, as defined as 
r =(xx + yy+ ZZ). (3.300) 


With the dot product evaluated in the exponential, the TM and TE 


waves are specified as 


EIM-p  gJko(xsinoszcoso) (x Coo zsing) (3.5 8) 


and 


¡A e JKo(x sina +ZCOSQ) y (3 5 9) 
s j 
The relative magnitude and phase of the computed field in the 
direction of the theoretical field are calculated and displayed. This 
"scaled dot product" is defined as 


Eee; 
E ' 





Ego = (3.60) 


where E is the computed field, e; is a unit vector pointing in the 


same direction as either the TM or TE incident field specified in 


Equations (3.58) and (3.59), and 


E; -E, e Jko& sina +zcosa) (3.61) 

The remaining field term that is computed and displayed by CORFLD 
is the error term. The error term is that portion of the computed 
field which is perpendicular to the incident field vector. In 
normalized form this becomes, 
E-(E»e;)e; 

E; 
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E = 


error 


(3.6 2) 





The relative magnitude of the error term is calculated and displayed 
by CORFLD as a measure of electromagnetic disturbance caused by 
the presence of the radome. This built-in theoretical field 
comparison is used in Chapter V as a measure of modelling accuracy 
by comparing the generated theoretical fields to the fields computed 
and assembled by EMRAD and CORFLD for an ideal radome 


constructed of air (e, = 1). 


C. CANTING THE PLANAR ARRAY 

The planar antenna within a radome can be canted in any 
direction allowed by the gimbal system that is used with the 
antenna. In order to evaluate the field components across a canted 
planar antenna within the radome core, the coordinates of each 
point on the antenna must be converted from the local planar 
coordinate system to the global radome core coordinate system. 
CORFLD uses a transformation from local planar to global core 
coordinates, assuming a two-axis gimbal system where the antenna 
is rotated about the original z-axis, as shown in Fig. 3.4, followed by a 
rotation about the new y-axis, as shown in Fig. 3.5. 

The transformation matrix is formed by multiplying two 
matrices, one for each axis rotation [Ref. 12]. The transformation for 


the rotation about the z-axis is provided by the matrix 


X 4cosp —sinB OO] |x’ 
y|=| sinB cosp 00| |y'|. (3.63) 
Z 00 00 10] Iz 
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Figure 3.4 Rotation about the z-axis 
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Figure 3.5 Rotation about the new y-axis 
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The transformation for the rotation about the y-axis is provided by 


the matrix 


x' COSY 00 siny X 
y |s9: 7010 1.0 00 Ve alae (3.64) 
z -siny OO cosy z" 


The complete transformation is given by the product 


X -cosp —sinB 0.0 COSY OO  siny} |x" 
y| =| sinB cosp 0.01 .| 00 1.0 00| |y" 
Z 0.0 0.0 1.0 -siny OO cosy} |z" 
cosycosD -sinB siny cosB| |x" 
= |cosysinB  cosB Siny sinB| |y"] . (3.00) 
-siny 00 COSY zi 


This transformation matrix is used by CORFLD to determine the 
global radome coordinates from the local antenna coordinates for a 
gimbaled antenna rotated by p degrees about the z-axis, then y 
degrees about the new y-axis. With this computational tool the fields 
that are received by an antenna within the core of the radome can be 
determine for any orientation of the antenna and for any computed 


incident field. 
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IV. PROGRAM EXECUTION AND DATA FLOW 


A. THE RADOME MODEL DESIGN PROCESS 

The design process developed by this thesis is represented by 
Fig. 4.1. The radome design procedure is an iterative process of 
selecting radome size and material parameters, radome shaping 
functions, and placement of the axis origin in order to allow for an 
optimal mesh formation by EMRAD. 

The flow of data starts when RADSHP is used to design the 
radome. RADSHP generates a material file which is ready for use by 
EMRAD and a structure data file which is in the format required by 
CURVE DIGITIZER. CURVE DIGITIZER is used to visually check the 
radome structure as designed in RADSHP. After verifying the 
structure data file, the program EMCADIN, developed by Connolly 
[Ref. 4], is used to translate the data file into the form required by 
EMRAD. After data translation, the program EMESH, developed by 
Morgan [Ref. 3], is used to check that the mesh to be generated in 
EMRAD is optimally placed. The placement of the origin in RADSHP 
when generating the structure data file can be modified by a 
placement factor. With the material and structure files are prepared, 
EMRAD is executed to generate the interior coefficients of expansion 
for use by CORFLD. A sample annotated execution through the data 


flow follows. 
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B. RADSHP EXECUTION 

RADSHP is a menu-driven program that accepts data from the 
keyboard. When the executable version of RADSHP is available, 
RADSHP is invoked by typing 


> RADSHP 


at the DOS prompt. The program will load itself, then provide 


introductory information. 


WELCOME TO RADSHP 


THIS PROGRAM GENERATES RADOME MATERIAL AND SHAPE DATA 
FILES FOR USE WITH EMRAD, BASED ON USER INPUT. 


RADSHP prompts for data file names for the material and structure 
output files. Sample keyboard inputs are provided. User inputs are 


preceded by "==>". 
TWO FILES ARE REQUIRED TO HOLD PROGRAM OUTPUT DATA. 
PLEASE ENTER THE NAME OF THE FILE TO BE USED TO 
HOLD THE PROGRAM MATERIAL PARAMETER DATA. THE EXTENSION 


MAT WILL BE ADDED AUTOMATICALLY, I. E. YOURNAME.MAT 


==> MATFIL 


PLEASE ENTER THE NAME OF THE FILE TO BE USED 
43 


TO HOLD THE PROGRAM OUTPUT DATA. THE EXTENSION .DAT WILL 
WILL BE ADDED AUTOMATICALLY, I. E. YOURNAME.DAT 


==> STRFIL 


RADSHP will allow the construction of radomes with up to five layers 
of construction. At present, each of these layers must completely 


enclose the previous layered region. 


ENTER THE NUMBER OF LAYERS (.LE. 5) 


=> 0 


RADSHP prompts for the parameters of each layer of construction. 
Each layer of the radome will have material constants, 
an outer height of the radome layer, H, 
a outer base radius, R, and 


a shape function, & shape factor parameters (if appl.), 


The first "layer" will be the inner core of the radome, 
normally filled with air. 


Layer No: 1 
This 1s the inner core of the radome 
Please enter the radome layer height, H, 


and the outer base radius, R, for this layer. 


==> 1.25,0.50 
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RADSHP provides a menu of the shaping functions available for 
radome construction. The curved parabola is not normally used for 
outer layers of the radome, but is included as an option for inner 


layer construction. 


The following shaping functions are available, 
1. General Ogive 

. Tangent Ogive 

. von Karman 

. Power Series 


. Parabola, curved 


OV tC 4d UC WN 


. Parabola, capped by pointed apex 


Indicate desired shaping function by entering a number 
1,2, ... or 6 


==> 2 


Each shaping function requires input of their respective 
shaping factors. These are prompted for and then listed for each 
layer as the radome is constructed. The tangent ogive and the von 
Karman functions have shaping factors that are determined by the 
height and base radius of the radome layer; these factors are 
determined and listed. The factor designation follows the form 


adopted in the equation in Chapter II. 


The tangent ogive shaping function determines two 
factors for the radome construction. These are calculated 
based on the height and radius of the radome layer. 
These factors are F4 = 6.50 and F5= 6.38 
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LAYER SHAPE HEIGHT RADIUS Fl F2 SES F4 ES F6 —B7 


] TAN Ogive 1.25 50 00  .00 .00 650 6.38 .00  .00 


Layer No: 2 


This is the outer layer of the radome 


Please enter the radome layer height, H, 
and the outer base radius, R, for this layer. 


==>" 130.035 


The following shaping functions are available, 
1. General Ogive 

. Tangent Ogive 

. von Karman 

. Power Series 


. Parabola, curved 


QV tA pw N 


. Parabola, capped by pointed apex 


Indicate desired shaping function by entering a number 
ZOO 


== 3 
The tangent ogive shaping function determines two 
factors for the radome construction. These are calculated 


based on the height and radius of the radome layer. 
These factors are F4 = 6.42 and F5 - 6.28 
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LAYER SHAPE HEIGHT RADIUS Fl FRPS F4 


F5 F6 F7 


- um un un Wn "VD VD "VID VD wm aa a "m m see es "un V -—— a on oo - am - a - m — -——— - oe eee = 


] TANOgive 1.25 .50 00 00 .00 6.50 6.38 .00  .00 


2 TAN Ogive 1.30 E35 00  .00 .00 6.42 6.28 .00 .00 


After all the structure forming data is entered, RADSHP 
prompts for the material relative permittivity parameter for each 


layer of construction. Layer 1 is the inner-most region and the 


permittivity will normally equal unity. 


Layer No: 1 


This is the inner core of the radome 
Enter Er for this Layer: 


—— NI 0 


Layer No: 2 


This is the outer layer of the radome 
Enter Er for this Layer: 


2-5040 
With the structure and material parameters entered, RADSHP 
prompts for headers to identify the radome constructed. This header 


information is stored in the material data file and will follow the data 


files throughout the execution of the programs in this package. 
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THE HEADERS ALLOW THE USER TO IDENTIFY THIS 
SET OF DATA FROM ALL OTHER SETS. 


PLEASE ENTER HEADER +1 (64 CHARACTERS MAX) 


==> This is a test data set for RADSHP example 


PLEASE ENTER HEADER +2 (64 CHARACTERS MAX) 


==> Header info stays with data set throughout 


With all the data input and calculations completed, RADSHP 
reports the generation of the material and structure files and 


terminates. 


Computations and output are completed, the data file, MATFIL.MAT 
holds the material specifications for input to EMCAD. 


The data file, STRFIL.DAT, holds the radome shape data for 
viewing in the CAD package, CURVE DIGITIZER. 


The structure file required for input to EMCAD can be 
generated by the executing the program EMCADIN, 
using, STRFIL.DAT, as input to EMCADIN, 


The output file from EMCADIN will be in proper format 
for input to EMCAD for the structure file. 


Stop - Program terminated. 
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C. RADOME CONSTRUCTION VERIFICATION. 

The construction of the radome by RADSHP is verified by 
entering the CURVE DIGITIZER package and displaying the structure 
data file generated by RADSHP. Figure 4.2 is a copy of the CURVE 
DIGITIZER display of the radome shape created in the example 
execution above. Upon review and acceptance of the structure data 
file, EMCADIN is used to translate the data file into the form required 
by EMRAD. EMCADIN is invoked by typing 


>EMCADIN 
at the system prompt. The program will load itself and displays the 


greeting heading and prompts for the structure data file as input. 
WELCOME TO EMCADIN 
THIS PROGRAM ACCEPTS PROPERLY FORMATTED INPUT DATA 
FROM CURVE DIGITIZER AND CONVERTS IT TO A FORM 
WHICH CAN BE USED BY EMCAD. 
PLEASE PRESS ENTER TO CONTINUE. 
INPUT FILENAME 
INPUTEN IS THE INPUT FILENAME OF THE FILE CONTAINING 
THE DATA OF INTEREST WHICH WILL BE INTERPOLATED ON AND 
CONVERTED TO A FORM WHICH CAN SUBSEQUENTLY BE USED BY 


EMCAD 


PLEASE ENTER THE NAME OF THE INPUT FILE. THE EXTENSION 
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Z-axis 


Xaax is 


Figure 4.2 Tangent ogive created by RADFLD 
as displayed by CURVE DIGITIZER 
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OF THE FILENAME MUST BE .DAT, LE. FILENAME.DAT 


==> STRFIL.DAT 


EMCADIN then prompts the user for the data file name to hold the 


Structure file that is prepared for EMRAD execution. 


OUTPUT FILENAME 


OUTPUTEN IS THE OUTPUT FILENAME OF THE FILE CONTAINING 
THE DATA WHICH HAS BEEN INTERPOLATED AND CONVERTED TO A 
FORM WHICH WILL BE USED BY EMCAD. 


PLEASE ENTER THE NAME OF THE OUTPUT FILE. THE EXTENSION 
OF THE FILENAME MUST BE INCLUDED I.E. FILENAME.DAT 


==> STRFIL.STR 


There are several data files used within these packages, the 
default file name extensions should be: ".MAT" for material data files, 
".DAT" for structure data files suitable for CURVE DIGITIZER, and 
"STR" for structure files suitable for EMRAD. 

After receiving the data file names to work with, EMCADIN 
prompts for angular resolution requirements. EMCADIN uses linear 
interpolation to provide the capability to generate output data with 
better resolution than the input data. The default resolution of data 
generation by RADSHP is one thousand points in the x-z plane, using 
180 points of angular resolution in EMCADIN is sufficient for good 


results. 


ul 


RESOLUTION 


DELTHE IS THE USER INPUT VALUE OF THE DESIRED THETA 
RESOLUTION IN DEGREES. 


PLEASE ENTER THE DESIRED DELTA THETA VALUE IN DEGREES. 


==> 180. 


Upon completion of the data file translation by EMCADIN, the 
data files required for execution by EMRAD are complete and in the 
proper format. MATFIL.MAT is the example material file and 
STRFIL.STR is the example structure file. Before executing EMRAD, 
the efficiency of mesh generation to be expected during EMRAD can 
be previewed by using the program EMESH developed by Morgan. 
Figure 4.3 shows the mesh generated by EMESH for the example 
radome, as displayed by CURVE DIGITIZER. 


D. EMRAD EXECUTION 

With the data files generated and an optimal mesh solution, the 
data flow is ready to initiate EMRAD. EMRAD is a menu-driven 
program developed by Connolly and Morgan and has been slightly 
modified to provide the interior spherical coefficients of expansion 


for use by CORFLD. EMRAD can be invoked by typing 


>EMCAD 


at the DOS prompt. The program will load itself and provide 
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Figure 4. 5 Tangential ogive with mesh generated by EMESH 
as displayed by CURVE DIGITIZER 
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introductory information. 


eee WELCOME TO EMRAD ****** 


THIS PROGRAM COMPUTES EM FIELD SCATTERING FROM PENETRABLE 
BODIES OF REVOLUTION USING THE COUPLED AZIMUTHAL POTENTIAL 
(CAP) FORMULATION IN CONJUNCTION WITH A VARIATIONAL FINITE- 
ELEMENT TECHNIQUE AND A TRI-REGIONAL UNIMOMENT METHOD. THE 
NUMERICAL SOLUTION IS PERFORMED USING A TWO-SWEEP RICCATI 
TRANSFORM WITH CIRCUMFERENTIAL MARCHING. 
UPDATED VERSION OF ORIGINAL U.C. BERKELEY CODE. 
MODIFICATIONS BY M.A. MORGAN MAR 1987 - DEC 1989 

E.M. CONNOLLY JUL 1987 - APR 1988 

R.J. VINCE SEP 1989 - DEC 1989 . 


PLEASE PRESS ANY KEY TO CONTINUE. 


ERA RR EMRAD INPUT A ale ok ok ale ale ale 


EMRAD ALLOWS A USER TO INPUT THE NECESSARY INFORMATION 
ACCORDING TO THE LEVEL OF HIS EXPERTISE. 


1 NOVICE LEVEL - GIVES BRIEF EXPLANATIONS/DESCRIPTIONS 
OF INPUT VALUES. DESCRIBES THE REQUIRED 
FORMAT FOR THE INPUT VALUE. GIVES 
TYPICAL VALUES WHERE APPLICABLE. 


2 EXPERT LEVEL - ASSUMES THE USER IS FAMILIAR WITH THE 
INPUT PARAMETERS AND FORMATS. SIMPLY 
PROMPTS FOR REQUIRED INPUTS. 


3 EXIT EMRAD 
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PLEASE SELECT THE LEVEL OF EXPERTISE BY ENTERING 1, 2, OR 3. 


The names of the data files generated by RADSHP and translated by 
EMCADIN are entered at the keyboard when prompted by EMRAD. 


EMRAD REQUIRES DATA INPUT FROM TWO FILES PREPARED 
IN ACCORDANCE WITH THE USERS MANUAL. THESE FILES MUST 
CONTAIN MATERIAL PARAMETERS AND STRUCTURE PARAMETERS, 
RESPECTIVELY. 


ENTER MATERIAL PARAMETER FILE (D:FILENAME.EXTENSION) : 


==> MATFIL.MAT 


ENTER STRUCTURE DATA FILE (D:FILENAME.EXTENSION): 


==> STRFIL.STR 


EMRAD opens both files and displays the header information which 
RADSHP placed in the material file. 


This is a test dataset for RADSHP example 
Header will stay with dataset throughout 


EMRAD creates several datafiles. The first will be the diary file for 
the execution of EMRAD. This file will have the ".OUT" extension. 


The second data file generated will have the same root name as the 


S5 


diary file, this is the data file which holds the spherical coefficients of 
expansion for the core field generation. This second file will have the 


extension ".INT" and will be used as input to CORFLD. 


THE OUTPUT DATA FILE IS THE FILE CONTAINING THE OUTPUT 
RESULTS OF ALL NUMERICAL CALCULATIONS CONDUCTED BY EMRAD. 
THE FORMAT FOR THIS INPUT IS FILENAME ONLY. NO EXTENSION IS 
REQUIRED OR DESIRED. EMRAD AUTOMATICALLY APPENDS AN 
EXTENSION OF .OUT TO YOUR FILENAME. 


A SECOND OUTPUT DATA FILE WILL BE CREATED TO STORE THE 
COEFFICIENTS FOR THE INTERIOR CORE EXPANSION. EMRAD WILL 
APPEND AN EXTENSION OF .INT FOR THIS DATA FILE. 

PLEASE ENTER THE FILENAME OF THE OUTPUT DATA FILES. 


==> INTFIL 


EMRAD was developed by Morgan and Connolly to determine the 
scattered field and the radar cross section of the body of revolution 
under study. Four output data files are created for input to separate 


graphing routines. 


THE OUTPUT GRAPHICS DATA FILE IS THE FILE CONTAINING THE 
OUTPUT DATA FROM EMRAD TO BE USED AS INPUT TO GRAPHING 
ROUTINES. THE FORMAT FOR THIS INPUT IS FILENAME ONLY. NO 
EXTENSION IS REQUIRED OR DESIRED. EMRAD AUTOMATICALLY 
APPENDS AN EXTENSION OF .TMT, .TMP, .TET, AND .TEP TO YOUR 
FILENAME AS IT PRODUCES FOUR OUTPUT FILES FOR GRAPHICS. 


.TMT ----- » TM INCIDENCE, F-THETA 
IMP ----- > TM INCIDENCE, F-PHI 
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.TET ----- » TE INCIDENCE, F-THETA 
TET ----- > TE INCIDENCE, F-PHI 


PLEASE ENTER THE FILE NAME OF THE OUTPUT DATA FILE. 


==> SCATGRF 


EMRAD prompts for a graphics caption to be included with the data 
set, this caption is passed along to CORFLD and is included with each 
plotted figure. This caption affords the user the ability to identify 


this data set from all other sets. 


THE GRAPHICS CAPTION IS A PERSONALIZED CAPTION ALLOWING THE 
USER TO IDENTIFY THIS SET OF GRAPHS FROM ALL OTHER SETS. 

THE MAXIMUM LENGTH OF THIS CAPTION IS 64 CHARACTERS. 

NOTE: WHEN USED WITH THE GRAPHICS PACKAGE, THE PROGRAM IS 
ABLE TO DIFFERENTIATE BETWEEN UPPER CASE AND LOWER CASE 
CHARACTERS. 


PLEASE ENTER ANY GRAPHICS CAPTION YOU DESIRE 


==> Graphics caption 


As discussed in Appendix A, EMRAD utilizes a semi-annular 
conformal mesh in order to conduct electromagnetic field 
calculations. EMRAD prompts the user for a minimum and maximum 
mesh density. EMRAD conducts the analysis of the radome structure 
using a linear variation of density from maximum density decreasing 


to minimum density. 


DMIN AND DMAX ARE PARAMETERS OF MESH DENSITY IN TERMS OF 
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ELEMENTS/INTERIOR LAMBDA. INPUT VALUES ARE EXPECTED TO BE 
REAL, I.E. THE DECIMAL POINT MUST BE INCLUDED. 
TYPICAL VALUES ARE DMIN = 10. AND DMAX = 15. 


PLEASE INPUT DMIN 


S 


PLEASE INPUT DMAX 


==> 18 


EMRAD allows for analysis of up to five incident fields. EMRAD 
prompts for the number of incident angles, and then for the value of 
each angle, in degrees. EMRAD uses the poynting vector direction for 
the incident angle. This translates to an angle of 180 degrees for an 
angle incident at the radome's appex, 135 degrees for an angle 
incident 45 degrees off the apex, and 90 degrees for an incident 


angle off the beam of the radome. 


THE NUMBER OF INCIDENT FIELD ANGLES IS THE TOTAL NUMBER OF 
INCIDENT FIELDS THAT IMPINGE ON THE OBJECT OF INTEREST. 

THIS PROGRAM ALLOWS A MAXIMUM OF FIVE INCIDENT FIELD 
ANGLES. WHEN ENTERING YOUR ANSWER PLEASE DO NOT INCLUDE A 
DECIMAL BECAUSE THE INPUT MUST BE IN INTEGER FORMAT, LE. 


NA =3 , NA=1 
PLEASE INPUT THE NUMBER OF INCIDENT FIELD ANGLES. 


zc 
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ENTER INC FLD ANGLE (DEG) FOR # 1 
==> 180 
ENTER INC FLD ANGLE (DEG) FOR # 2 


==> 135 


EMRAD performs scattering calculations by starting at @ = O, and 
progressing clockwise to @ = 180 in equal increments, the menu 


allows the user to determine the size of the increments. 
THE NUMBER OF SCATTERING FIELD THETA POINTS DETERMINES THE 
SPACING BETWEEN THETA POINTS DURING EMRAD ITERATIONS AND 
CALCULATIONS. 
DELTA THETA = 180/ (NUMBER THETA POINTS - 1) SO... 
NUMBER THETA POINTS = 37 ----- > DELTA THETA = 5 DEGREES 
NUMBER THETA POINTS = 19 ----- > DELTA THETA = 10 DEGREES 
WHEN ENTERING YOUR ANSWER, PLEASE DO NOT INCLUDE A DECIMAL 
BECAUSE THE INPUT MUST BE IN INTEGER FORMAT. I.E. NT = 19 


PLEASE INPUT THE NUMBER OF SCATTERING FIELD THETA POINTS 


==> 19 


As part of the scattering analysis, EMRAD allows for eight angles of 
user observation. 
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THE NUMBER OF PHI ANGLES IS THE TOTAL NUMBER OF PHI ANGLES. 
THIS PROGRAM ALLOWS A MAXIMUM OF THREE PHI ANGLES. WHEN 
ENTERING YOUR ANSWER PLEASE DO NOT INCLUDE A DECIMAL 
BECAUSE THE INPUT MUST BE IN INTEGER FORMAT. LE. 


NP=3 , NP=1 


PLEASE INPUT THE NUMBER OF PHI ANGLES. 


EMRAD truncates the modal analysis of the electromagnetic field 
infinite fourier expansions for the incident plane waves and the 
interior core as well as the exterior scattered fields. Mathematical 
estimates of truncations are provided and the user is prompted for 
truncation limits. Higher truncations may increase the modelling 


accuracy. 


ENTER MSTOP (.LE. 13) 
ESTIMATED "MINIMUM" VALUEIS: 4 


==> 5 


Modal truncations for the number of spherical harmonic modes in 
the spherical radome core and external scattering region are 
estimated. The user is prompted for truncation limits for the 


spherical harmonics. 
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ENTERING NO. OF INTERNAL AND EXTERNAL EXPANSION MODES 
ESTIMATED (KI*RMIN) "MINIMUM" INTERNAL N1= 4 
ENTER N1 (GE. MSTOP): 


==>. 5 


ESTIMATED (KO*RMIN) "MINIMUM" EXTERNAL N2= 9 
ENTER N2 (.GE. MSTOP): 

===> 9 
Check that N1+N72 is .LE. 35 .... Otherwise Abort 


==> [RETURN] 


When executing the numerical solution, EMRAD requires a large 
amount of temporary disk storage space. The amount of storage 
required is estimated and provided to the user. The user must select 
a disk drive that has adequate memory available for usage. The use 
of a ram drive configuration can increase the speed of execution by 


EMRAD. 


***** Temporary Storage of Riccati Matrices ***** 
Estimated Disk Space Needed in MegaBytes: 1.930000E+00 
ENTER Hard Drive or RAM Disk Letter (C,D, E... ): 


E UE 
6 1 


With all of the data and parameters specified, EMRAD begins the 
numerical analysis of the electromagnetic interaction of the radome 
Structure and the specified incident plane waves. Execution is 


monitored by reports to the terminal screen. 
INDEXING PROGRAM THROUGH VALUES OF M 


M-LOOP ..M= 0 

CALL MESH 

ENTER I-LOOP; NO STEPS TO COMPLETE: 38 
I= 1  NOSTEPS TO GO: 37 

CALL LODER 

CALL MARCH 

I= 2 NO STEPS TOGO: 36 


EX M-LOOP, DATOUT 
aaao k EMR AD COMPLETED eo %44k kkk kk kkk 


Stop - Program terminated. 


Upon termination of the program EMRAD, the coefficients for analysis 
of the radome core fields have been stored in the data file, ready for 


CORFLD execution. 


E. CORFLD EXECUTION 
CORFLD is a menu driven program that accepts data from the 


data file created by EMRAD and from the keyboard. When the 
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executable version of CORFLD is available, CORFLD is invoked by 
typing 


> CORFLD 


at the DOS prompt. The program will load itself, then provide 


introductory information. 


WELCOME TO CORFLD 


CORFLD ASSEMBLES THE FIELD ON THE INTERIOR OF A RADOME 
WITH THE COEFFICIENTS SUPPLIED BY AN EMRAD OUTPUT FILE. 


CORFLD prompts for the name of the data file generated by EMRAD. 
Only the root file name is required, CORFLD attaches the ".INT" to the 


root file name. 


THE OUTPUT DATA FILE IS THE FILE CONTAINING THE OUTPUT 
RESULTS OF INTERIOR FIELD CALCULATIONS CONDUCTED BY EMRAD. 
THE FORMAT FOR THIS INPUT IS FILENAME ONLY. 

NO EXTENSION IS REQUIRED OR DESIRED. 

CORFLD AUTOMATICALLY APPENDS AN EXTENSION OF .INT TO YOUR 
FILENAME. 


PLEASE ENTER THE FILENAME OF THE OUTPUT DATA FILE. 


==> INTEIL 


Header information that has followed the data flow is obtained from 


the ".INT" data file by CORFLD and displayed on the screen. Notice 
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that the graphics caption entered during EMRAD execution as well as 


the headers entered during RADSHP execution are displayed. 


TEST DATA FOR USER EXAMPLE 
THIS HEADING FOLLOWS THE DATA THROUGHOUT 
Graphics caption 


CORFLD reads the material parameters stored in the ".INT" file and 


displays them on the screen. 


SRE- (1.000000,-5.000000E-07) 
SRU- (1.000000,-5.000000E-07) 
KR = (1.000000,-1.000000E-06) 


CORFLD reads the number of incident angles evaluated by EMRAD 
and displays each angle as the poynting vector angle of incidence and 
also as the angle of incidence to the radome exterior. The radome 
incident angle is displayed so that the user can use the angle as a 


reference when canting the planar antenna. 


Incoming Inc. Angle + l 
EMRAD Poynting angle of incidence = 180.000000 
Radome incident angle = 0.000000E+00 


Incoming Inc. Angle # 2 
EMRAD Poynting angle of incidence= 135.000000 
Radome incident angle = 45.000000 


After displaying the angles of incidence available for consideration, 


EMRAD prompts the user for the angle of interest. Only one angle 
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may be selected for each execution of CORFLD. In order to consider 


other angles, CORFLD must be re-executed by the user. 


Select Angle of Interest By Number 


The angle of consideration is angle + 1 =  0.000000E+00 degrees 


After an angle is selected for analysis, the user may enable a rotation 


of the planar array. 


The planar array is set by default to lie in the 
x-y plane, with boresight in the z-direction. 


Do you want to change the array boresight ? 
( e Or ENP ) 


==> N 
CORFLD reads the data from the ".INT" file and assembles the interior 


field components. 


MSTART= O MSTOP= 5 
Reading Coefficients and Constructing 
Fields for Each Azimuthal "m" Mode 


Terminating M-Loop at m=1 for +/- Z-Axis Incidence 


Azimuthal Mode "m" = 0 
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Reading Coefficients 
Azimuthal Mode "m" = ] 
Reading Coefficients 


Generate Fields Across Array 
SPHERICAL FIELDS COMPUTED 


When the interior fields have been computed, the user may check for 


the field component values for any point on the planar array. 


Check Values of Coordinates and E-Field ? (Y/N): 


mo 

Enter Array Triple Index to Check (J,L,D 
J= x-point 1,2,.., 50 
L= y-point 1,2,.., 50 


I= 1 for TM, 2 for TE 


XP,YP,ZP: -5.115592E-02 -5.115592E-02  0.000000E-00 
XC,YC,ZC: -5.115592E-02 -5.115592E-02 0.000000E+00 


IE-RADI 2 5.023160E-01 
IE-THI 2 1.391270E-02 
IE-PHII 2 5.022836E-01 
IE-TOT! - 7.104954E-01 


Look at Another Point ? (Y/N): 


== Ww 


The user is supplied with the selection menu to allow access to any 


component representation of the interior fields across the planar 
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antenna which is oriented as specified by the user. Component phase 
and magnitudes of the interior fields as well as the total field 


magnitudes are available, 


CORFLD OUTPUT GRAPHICS PRESENTATION 
CORE FIELD REPRESENTATION SELECTION MENU 


Please select the type of field representation of interest 
CORE will display the following field representations: 


TM INCIDENT FIELD, E-THETA 

TM INCIDENT FIELD, E-PHI 

TM INCIDENT FIELD, E-RADIAL 

TM INCIDENT FIELD, E-TOTAL FIELD MAGNITUDE 


E 


5. TE INCIDENT FIELD, E-THETA 

6. TE INCIDENT FIELD, E-PHI 

7. TE INCIDENT FIELD, E-RADIAL 

8. TE INCIDENT FIELD, E-TOTAL FIELD MAGNITUDE 
9. COMPARE COMPUTED FIELD TO INCIDENT FIELD 
10. CHANGE ASPECT RATIO 

0. FINISHED WITH THIS ANGLE OF INCIDENCE 


Indicate your selection entering "1","2",...,"10",or "O" 


4 


When the user has made a selection from the menu, the three 


dimensional (3-D) plotting routine is called and the user is allowed to 
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adjust the view angle for best observation of the 3-D plot, see 


Fig. 4.4, 
CALL PLOTTING ROUTINE 


A black and white 'stick' representation as well as a color 3-D fill 
representation are available to the user. The ‘stick’ representation 


generates the clearest output for printing and is used for all the data 


presented here. 


Enter 0 for 3D-stick or 1 for 3D-fill option 


(** SEE FIGURE 4.4 **) 


Default view angles are: phi=-45.,theta=70. 
ENTER NEW ANGLES PHI & THETA (DEG) OR (9,9) WHEN DONE 


==> 9,9 


When the user is finished with an electric field component, the new 


view angle "9,9" is entered and the user is returned to the core field 


selection menu. 


CORE FIELD REPRESENTATION SELECTION MENU 


9. COMPARE COMPUTED FIELD TO INCIDENT FIELD 
Indicate your selection entering "1","2",...,"10",or "O" 
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The datafile used to generate fields is intfil. INT 
TEST DATA FOR USER EXAMPLE 
THIS HEADING FOLLOHS THE DATA THROUGHOUT 
Graphics caption 
TMH INCIDENT FIELD, E-TOTAL FIELD MAGNITUDE 
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min. Z= +.631; max. Z = *.825; average Z - +.656 
Inc. Á = +,888;: Plane G = +.88B; View: phi=+45.8 theta=+75.80 


Figure 4.4 3-D plot of E-total field magnitude for example ogive radome 
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When the user wishes to see the amount of perturbation caused to 


the incident plane wave by the radome, selection "9" is entered, 


~ 


Selection "9" sends the user into Menu 2, for core field comparison to 
the ideal theoretical plane wave which would be incident on the 


antenna without the radome present. 


CORE FIELD REPRESENTATION SELECTION MENU 2 


Please select the type of field representation of interest 
CORE will display the following field representations: 


TM INCIDENT FIELD, Computed Field dot Ideal Field 
TM INCIDENT FIELD, Scaled Error component 
TE INCIDENT FIELD, Computed Field dot Ideal Field 
TE INCIDENT FIELD, Scaled Error component 


PWN — 


Indicate your selection entering "1","2","3", or "4" 
—— 2; 


CALL PLOTTING ROUTINE 
Enter O for 3D-stick or 1 for 3D-fill option 


Default view angles are: phi=-45.,theta=70. 


(** SEE FIGURE 4.5 **) 
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The datafile used to generate fields is intfil. INT 
TEST DATA FOR USER EXAMPLE 
THIS HEADING FOLLOHS THE DATA THROUGHOUT 
Graphics caption 
TM INCIDENT FIELD, Scaled Error component 
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= 4.013; max. Z = +.436; average Z = *.221 
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Figure 4.5 3-D plot of scaled error component for example ogive radome 
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ENTER NEW ANGLES PHI 8 THETA (DEG) OR (9,9) WHEN DONE 
==> 90 


CORFLD OUTPUT GRAPHICS PRESENTATION 
CORE FIELD REPRESENTATION SELECTION MENU 


0. FINISHED WITH THIS ANGLE OF INCIDENCE 
Indicate your selection entering "1","2",...,"10",or "O" 
==>() 


Stop - Program terminated. 


When the user is finished with the angle of incidence selected, "O" is 
selected at the selection menu and the program execution is 


terminated. 


te 


V. PROGRAM VALIDATION AND RESULTS 


In order to validate the execution of EMRAD and CORFLD, the 
tangent ogive radome generated in Chapter IV was modified to have 
a relative permittivity equal to that of air, e,=1.0. The field 
coefficients were generated by EMRAD and the core fields were 
computed by CORFLD and compared to the unperturbed incident field 
as listed in Chapter III. 

Figures 5.1 and 5.2 show the magnitude and phase of the 
scaled dot product, as defined in Chapter III, between the computed 
and theoretical values for the TM and TE incident fields, respectively, 
for a "nose-on" incident angle. Each comparison shows a 94.9% 
magnitude matching with an almost constant 2.59° phase differential 
across the array. Figure 5.3 exhibits the scaled error component for 
the TM and TE incident fields, with an average magnitude of error 
component of 0.4%. The comparisons for this "nose-on" angle of 
incidence show that there is good agreement between the theoretical 
and calculated field components. It is observed that the TM and TE 
components displayed are 90% out of phase with each other. This is 
expected for "nose-on" incidence since the TM wave has its direction 
component in the -x direction, while the TE  wave's direction 
component is in the y direction. 

Figure 5.4 shows the magnitude and phase of the scaled dot 


product between the computed and theoretical values for the TM 
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The datafile used to generate fields is airfil.1NI 
Shaping function = TAN ogive 
Material of formation is air, Er = 1 


TH INCIDENT FIELD, Computed Field dot Ideal Field, MAGN. 





min. Z = +.982; max. Z = +.996; average Z = +.949 
Inc. Á = +,880;Plane G = +.088; View: phi=+45.8 theta=+75.8 


Ihe datafile used to generate fields is airfil. 1NIT 
Shaping function = TAN ogive 
Material of formation is air, Er = 1 


TH INCIDENT FIELD, Computed Field dot Ideal Field, PHASE 
9. 08 





min. Z = +1.73; max. Z = +2.95; average Z = +2.59 
Inc. A = +. 888; Plane G = +.880; View: phi=+88.0 theta=+75.8 


Figure 5.1 Dot product magnitude and phase for TM 
computed fields determined for radome composed of air 
compared with theoretical plane waves 
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The datafile used to generate fields is airfil. INT 
Shaping function = TAN ogive 
Material of formation is air, Er = 4 


TE INCIDENT FIELD, Computed Field dot Ideal Field, MAGN. 
1.56 





min. Z = +.982; max. Z = +.996; average Z = +.949 
Inc. A = +. 800: Plane G = +.000; View: phi=+45.8 theta=+75.8 


The datafile used to generate fields is airfil. INT 
Shaping function = TAN ogive 
Material of formation is air, Er = 1 


TE INCIDENT FIELD, Computed Field dot ldeal Field, PHASE 
3. 00 


4. 88 
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Figure 5.2 Dot product magnitude and phase for TE 
computed fields determined for radome composed of air 
compared with theoretical plane waves 
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The datafile used to generate fields is airfil. 1NT 
Shaping function = TAN ogive 
Material of formation is air, Er = 1 


TE INCIDENT FIELD, Scaled Error component 
0.18 





min. Z = +. 000; max. Z = *.023; average Z = +. 004 
lnc. A = +. 8000; Plane G = +. 000; View: phi=+10.0 theta=+80. 8 


The datafile used to generate fields is airfil. INT 
Shaping function = TAN ogive 
Material of formation is air, Er = 1 


TM INCIDENT FIELD, Scaled Error component 
8.18 
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min. Z = +. 808080; max. Z = *.023; average Z = +.BB4 
Inc. A = +.080:Plane G = +.880; View: phi=+75.B theta=+88.8 


Figure 5.3 Scaled error component for TM and TE 
computed fields determined for radome composed of air 
compared with theoretical plane waves 
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The datafile used to generate fields is airfil. INT 
Shaping function = TAN ogive 
Haterial of formation is air, Er = 1 


TH INCIDENT FIELD, Computed Field dot Ideal Field, MAGN. 
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min. Z = +.540; max. z- +1.66; average Z = +.874 
Inc. A = +68.8; Plane G = +68.8; View: phi=+45.8 theta=+75.0 





The datafile used to generate fields is airfil. INT 
Shaping function = TAN ogive 
Material of formation is air, Er = 1 


TM INCIDENT FIELD, Computed Field dot Ideal Field, PHASE 
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min. Z = -27.9; max. Z = *60.1; average Z - -2.84 
Inc. A = *68.8: Plane G = +68.8; View: phi=+45.89 theta=+75.80 


Figure 5.4 Dot product magnitude and phase for TM computed fields 
determined for radome composed of air compared with theoretical 
plane waves for incident angle of 60 degrees 
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incident field, for an incident angle 60° from the "nose-on” aspect. 
The magnitude of the scaled dot product between the computed and 
theoretical TM incident wave shows modal variation around the 
antenna plane. The magnitude of the computed value in the 
direction of the theoretical plane wave varies between 54.0% and 
166%, with an average magnitude of 87.4%. The phase component of 
the dot product varies across the plane from -27.99 to 460.19, with 
an average phase of -2.849. The maximum phase differential occurs 
at the edge of the antenna plane closest to the point of wave 
incidence, as the incident wave is defined as arriving in the -x - z 
axis. This variation seen in the magnitude and phase of the scaled 
dot product in the direction of the theoretical, undisturbed plane 
wave is confirmed in the scaled error component of Fig. 5.5 with an 
amplitude variation from 8.4% to 88.6%, averaging at 23.6% of the 
computed field evaluated as an error term perpendicular to the 
theoretical fields expected. 

Figures 5.6 and 5.7 exhibit the behavior of the calculated TE 
60° incident plane wave, which is similar, but less in magnitude 
compared to the calculated TM wave. Less modal variation is seen 
along the exhibited plane of the antenna. The scaled error 
component has a symmetry about the x-axis. The scaled error 
amplitude varies form 0.3% to 33.1%, with an average of 12.5%. 

The form of the error components as analyzed for this radome 
constructed of air identifies a phase error in the determination of the 
computed fields. The data computations of the EMRAD subroutines 


have been previously validated by Morgan. These validations 
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The datafile used to generate fields is airfil. INT 
Shaping function = TÁN ogive 
Material of formation is air, Er = 1 


TM INCIDENT FIELD, Scaled Error component 
. 98 
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min. Z = +.B84;: max. T = Mc: average Z - 1.236 
Inc. A = +68.B;Plane G = +68.8; View! phi=+45.8 theta=+75.0 


Figure 5.5 Scaled error component for TM computed fields 
determined for radome composed of air compared with 
theoretical plane waves for incident angle of GO degrees 
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The datafile used to generate fields is airfil. INT 
Shaping function = TÁN ogive 
Material of formation is air, Er = 1 


TE INCIDENT FIELD, Computed Field dot Ideal Field, MAGN. 





min. Z - *.860; max. Z = +1.34; average Z = +.981 
Inc. A = +60.8; Plane G = +60.0; View: phi=+45.0 theta=+75.8 


The datafile used to generate fields is airfil. INT 
Shaping function = TAN ogive 
Material of formation is air, Er = 1 


TE INCIDENT FIELD, Computed Field dot Ideal Field, PHASE 
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min. Z - -18.2; max. Z = *4.96; average Z - -.43U 
Inc. f - *60.0; Plane G - *60.0; View: phi=+45.8 theta=+75.08 


Figure 5.6 Dot product magnitude and phase for TE computed fields 
deter mined for radome composed of air compared with theoretical 
plane waves for incident angle of 60 degrees 
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The datafiie used to generate fieids is airfii. INT 
Shaping function = TAN ogive 
Hateriai of formation is air, 


Er = 1 


TE INCiDENT FiELD, Scaied Error component 
6.56 
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+68.8; View: 
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Figure 5.7 Scaled error component for TM computed fields 
determined for radome composed of air compared with 


theoretical plane waves for incident angle of 60 degrees 
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considered the magnitude and phase of the scattered fields [Refs. 3 
and 4], and the magnitude of the core fields [Ref. 10] along the main 
coordinate axes. This work is the first to plot both the phase and 
magnitude variations across the entire two-dimensional plane and 
identifies the possible need for increased computational accuracy in 
the numerical routines of EMRAD. There is also the possibility of 
unresolved programming errors in CORFLD which should be further 


investigated. 


B. COMPUTATION RESULTS FOR VARIOUS RADOME MATERIALS 
EMRAD and CORFLD executions for the tangent ogive standard 
established in Chapter IV were run for varying material permittivity 
according to the material properties listed in Table 2.1. The 
comparisons of the computed fields to the unperturbed theoretical 
incident plane waves are included only for the "nose-on" incident 


angle which has good results for the ideal "air" construction. 


Fiber-reinforced glass, with e€e,=4.0; glass ceramic, with e,=5.5; and 
refractory oxide, with e€,=8.0, are included. The form of field 


comparisons are similar in shape for each radome considered, only 
the magnitude of the plots varies. As expected from the discussion 
in Chapter II, the attenuation of the incident plane wave, received at 
the plane of the antenna in the core, increases with increasing 
relative permittivity. When the radome was constructed of air, the 
computed magnitude of each incident wave was 94.9%. When 


fiber-reinforced glass construction is considered, Figures 5.8, 5.9, and 
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5.10, the average magnitude of the computed wave is 60.1% and the 
average scaled error component is 22.7%. For the radome 
constructed of glass ceramic, Figures 5.11, 5.12, and 5.13, the average 
magnitude is 66.9% and the average scaled error component is 13.3%. 
For the refractory oxide radome, Figures 5.14, 5.15, and 5.16, the 
average magnitude of received field is 87.6% and the average scaled 


error component is 20.6%. 
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The datafile used to generate fields is P4FIL. 1NT 
Shaping function = TAN ogive 
Material of formation is fiber-reinforced glass 


TH INCIDENT FIELD, Computed Field dot Ideal Field, MAGN. 
1.58 
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The datafile used to generate fields is p4fil. INT 
Shaping function = TAN ogive 
Material of formation is fiber-reinforced glass 


TH INCIDENT FIELD, Computed Field dot Ideal Field, PHASE 
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Figure 5.8 Dot product magnitude and phase for TM computed fields 
deter mined for radome composed of fiber-reinforced glass compared 
with theoretical plane waves for incident angle of 0 degrees 
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The datafile used to generate fields is p4fil. INT 
Shaping function = TAN ogive 


Material of formation is fiber-reinforced glass 


TE INCIDENT FIELD, Computed Field dot Ideal Field, MAGN. 
1.50 


1.23 
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The datafile used to generate fields is p4fil. INT 
Shaping function = TAN ogive 
Material of formation is fiber-reinforced glass 


TE INCIDENT FIELD, Computed Field dot Ideal Field, PHASE 





min. Z = -71.1: max. Z = -27.4: average Z = -51.8 
Inc. A = +. 800; Plane G = +.9808; View: phi=+88.8 theta=+88. 0 


Figure 5.9 Dot product magnitude and phase for TE computed fields 
determined for radome composed of fiber-reinforced glass compared 
with theoretical plane waves for incident angle of 0 degrees 
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The datafile used to generate fields is p4fil. INT 
Shaping function = TAN ogive 
Material of formation is fiber-reinforced glass 


TH INCIDENT FIELD, Scaled Error component 
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The datafile used to generate fields is p4fil. INT 
Shaping function = TAN ogive 
Material of formation is fiber-reinforced glass 


TE INCIDENT FIELD, Scaled Error component 





min. Z = +.813; max. Z = +.436; average Z = +.227 
Inc. A = +. 8000; Plane G = +. 8000; View: phi=+10.0 theta=+70.0 


Figure 5.10 Scaled error component for TM and TE computed fields 


determined for radome composed of fiber-reinforced glass compared 
with theoretical plane waves for incident angle of 0 degrees 
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The datafile used to generate fields is c55f il. INT 
TEST DATA FOR THESIS 


Radome is tangential ogive, 
constructed of ceramic, Er = 5.5 


TH INCIDENT FIELD, Computed Field dot Ideal Field, MAGN. 
1.50 


1.25 


atte 9. CS $3 
* SN IE 17» ae 
iu AR. 





min. Z = +,668:;: max. Z = +.883: average Z = +.669 
Inc. & - *.808008:;Plane G - *,.888: View: phiz*45.0 theta-*75.0 


The datafile used to generate fields is cSSfil.INT 
TEST DATA FOR THESIS 
Radome is tangential ogive, 


constructed of ceramic, Er = 5.5 
TIM INCIDENT FIELD, Computed Field dot Ideal Field, PHASE 
28.00 
8.58 8.88 -0.58 -8.50 8.88 8.58 
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Figure 5.11 Dot product magnitude and phase for TM computed fields 
determined for radome composed of glass ceramic compared 
with theoretical plane waves for incident angle of 0 degrees 
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The datafile used to generate fields is c5Sf il. INT 
TEST DATA FOR THESIS 
Radome is tangential ogive, 
constructed of ceramic, Er = 5.5 
TE INCIDENT FIELD, Computed Field dot Ideal Field, MAGN. 


1.58 
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min. Z = +.668; max. Z E *.803; average Z = +.669 
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The datafile used to generate fields is c55fil.INT 
TEST DATA FOR THESIS 
Radome is tangential ogive, 
constructed of ceramic, Er - 5.5 
TE INCIDENT FIELD, Computed Field dot Ideal Field, PHASE 


28.80 
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Inc. A = +.888;Plane G = +.888; View: phi=+45.8 theta=+75.8 


Figure 5.12 Dot product magnitude and phase for TE computed fields 


determined for radome composed of glass ceramic compared 
with theoretical plane waves for incident angle of 0 degrees 
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The datafile used to generate fields is c55f il. INT 
TEST DATA FOR THESIS 
Radome is tangential ogive, 
constructed of ceramic, Er = 5.5 
IM INCIDENT FIELD, Scaled Error component 
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The datafile used to generate fields is c55fil. INT 
TEST DATA FOR THESIS 
Radome is tangential ogive, 
constructed of ceramic, Er = 5.5 
TE INCIDENT FIELD, Scaled Error component 





min. Z = *.005; max. Z = +.312; average Z = +.133 
Inc. A = +.66@; Plane G = +.000; View: phi=#18.0 theta=+70.8 


Figure 5.13 Scaled error component for TM and TE computed fields 
determined for radome composed of glass ceramic compared 
With theoretical plane waves for incident angle of 0 degrees 
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The datafile used to generate fields is OX8FIL. INT 
Shaping function = TAN ogive 
Material of formation is refractory oxide, Er = 8 
TM INCIDENT FIELD, Computed Field dot Ideal Field, MAGN. 
2.008 





Z = +1.88; average Z = +.876 


min. 2 = +.634; max. 
Inc. Á = +.000:Plane G = +.000; Uieu: phiz*45.0 theta-*75.0 


The datafile used to generate fields is OX8FIL.INT 


Shaping function = TAN ogive 
Material of formation is refractory oxide, Er = 8 


IM INCIDENT FIELD, Computed Field dot Ideal Field, PHASE 
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Figure 5.14 Dot product magnitude and phase for TM computed fields 
determined for radome composed of refractory oxide compared 
with theoretical plane waves for incident angle of 0 degrees 
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The datafile used to generate fields is OX8FIL., INT 


Shaping function = TAN ogive 
Material of formation is refractory oxide, Er = 8 


TE INCIDENT FIELD, Computed Field dot Ideal Field, MAGN. 
2.868 
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The datafile used to generate fields is OXSFIL.INT 


Shaping function = TAN ogive 
Material of formation is refractory oxide, Er = 8 


TE INCIDENT FIELD, Computed Field dot Ideal Field, PHASE 
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Figure 5.15 Dot product magnitude and phase for TE computed fields 
determined for radome composed of refractory oxide compared 
with theoretical plane waves for incident angle of O degrees 
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The datafile used to generate fields is OX8FIL. INT 
Shaping function = TAN ogive 
Material of formation is refractory oxide, Er = 8 


™ INCIDENT FIELD, Scaled Error component 
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The datafile used to generate fields is OX8FIL.INT 
Shaping function = TAN ogive 
Material of formation is refractory oxide, Er = 8 


TE INCIDENT FIELD, Scaled Error component 
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Figure 5.16 Scaled error component for TM and TE computed fields 
determined for radome composed of refractory oxide compared 
with theoretical plane waves for incident angle of 0 degrees 
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VI. CONCLUSIONS 


This thesis has presented a useful set of programs to 
numerically model a radome, designed of specified shape, 
construction, and material, and to evaluate the perturbation of the 
electromagnetic field as it passes through the modelled radome to 
the antenna within the radome structure. The radome design 
methodology was discussed in Chapter Il, along with the radome 
shaping functions, wall construction, and choice of materials for 
radome construction. The program RADSHP is introduced and the 
shaping function, as well as the shaping factors for each function are 
explained. Chapter III developed the electromagnetic theory 
required to assemble the field components within the interior 
radome core, using the coefficients provided by the EMRAD program. 
The theory required for transforming the coordinates of the canted 
antenna plane to the core coordinates for an assumed two-axis 
gimbal was also discussed. Chapter IV explained the data flow for 
the radome model design process by using an example of the 
program executions for the design of a monolithic tangent ogive, with 
computational analysis by EMRAD and field component generation by 
CORFLD. The verification of the radome design and optimal placing of 
the core origin by the use of CURVE DIGITIZER software is also 
discussed. Chapter V discussed the computational errors discovered 


in the 3-D plotting of the field components. The scaled magnitude of 
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the computed perturbed plane waves for the "nose-on" incident 
angle agreed within 95% of the theoretical value and the phase of the 
computed field was an almost constant 2.59% out of phase with the 
theoretical incident plane wave. The errors increase when the angle 
of incidence is off the "nose-on" aspect, with modal amplitude 
variations from 54% to 166% of the theoretical plane wave. Further 
diagnostics are required to determine the source of this error. 

This is only the first stage in the development of a 
mathematical model for the interaction between the radome and the 
antenna system that the radome is designed to enclose. The model, 
in its present form, does not consider the presence of the antenna 
Structure within the radome core. Further work is needed to 
determine the effect of the antenna structure and its interaction with 
the radome. When such a model is developed the predicted data can 
be compared to data collected by a physical antenna system. The 
culmination of this effort could be the analytic development of 
corrective lenses, as referred to in Chapter II, to correct or minimize 
the perturbation of the incident electromagnetic field caused by the 


radome structure. 
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Appendix A 


EMRAD FIELD THEORY 

The program EMRAD uses the unimoment method mentioned in 
Chapter II, as developed by Morgan and Mei [Ref. 5] to decouple the 
interior and exterior spherical solutions. The expansion coefficients 
for the electromagnetic fields within the radome structure are 
generated as part of the scattering solution. CORFLD uses these 
coefficients for construction of the interior fields. The theory 
required by CORFLD to formulate the electromagnetic fields within 
the radome is discussed in Chapter II. 

The code developed by Morgan for EMRAD uses an optimized 
variational finite-element algorithm, in conjunction with a 
tri-regional unimoment method, to provide scattering solutions for 
each of multiple incident fields impinging on the specified structure. 
EMRAD represent a significant improvement in computational speed 
and versatility over computer codes which find the solution to the 
electromagnetic scattering and penetration problem by use of surface 
current integral equation formulations. This advantage results 
because of the highly sparse matrices embodied in the finite element 
solution vis-a-vis full matrices found using an integral equation 
formulation. An overview of the finite element scattering algorithm, 
along with the concept of coupled azimuthal potentials and the 
unimoment method are included in this appendix. Further 


information is available in references 3, 4, 5, 6, and 7. 
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A. COUPLED AZIMUTHAL POTENTIALS 
In the radome construction the material constitutive 
parameters, e(r) and u(r), are invariant to the azimuthal coordinate 4 
due to the axial symmetry of the radome. This axisymmetry allows 
implementation of vector field generation via coupled azimuthal 
potentials, as developed by Morgan, Chang and Mei [Ref. 7]. 
For ease of solution and representation, a normalized 
cylindrical coordinate system is adopted, as defined by (R , Z , 0) = 


(k, p, k, z, 0) where (p, z, 6) are standard circular cylindrical 
coordinates and k, = 2n/A, is the free-space wavenumber of the 


time-harmonic field. The equations in this form are related to the 
standard spherical coordinates (r, 8,6) as shown in Fig. A.l using the 


substitutions 


R = kọ r sin (A.1) 
and 
L= ki T coson. (A.2) 


The total electromagnetic fields are decomposed into azimuthal 


modes by an exponential Fourier series in the $-coordinate such that 


eoo 


E(RZ4) - Y, eu (R2) e"? (A.3) 
and 
n HZ - Y, h,(RZ e , (A.4) 
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Figure A.1 Rectangular (X, y, 2), spherical (r 8 4), and 


normalized cylindrical (R, Z,8) coordinate systems 
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where n, - 120 n Q is used to normalize the magnetic fields, so that 
h = no H. 

Using the axial symmetry of the radome construction, the 
modal fields decouple when the field expansions are substituted into 
Maxwell's equations, yielding [Ref. 7], after judicious cross 
substitutions the following modal field generating equations in terms 


of two scalar potentials, P,( R,Z,m ) and P,( R,Z,m ): 


do xen =j fn ( moxVE, — Rp, VI), ) (A.5) 
| — 

9 Cm S R (A.6) 

oxh, = jfm( moxVI, + Re, VI, ) (A.7) 
| -b 

p- hm = 2 (A.8) 


where the two dimensional gradient operator is defined as 


V -[R (80/or) * Z (83/02) ] (A.9) 


and f, is a multiplicative function given by 


| 
f (R2 = | 1 (RZ) € (RZ R=m? |. (A.1 0) 
Note that the potentials, [, and I, are proportional to the azimuthal 


components of the electric and magnetic fields, e,, and h,,. These 


m 
potential functions are continuous everywhere, including at dielectric 
and magnetic interfaces. This property of uniform continuity in the 


field is used to obtain a numerical solution of the fields. 


98 





The total fields may be described in term of the uniformly 


continuous scalar potentials T,( R,Z,m ) and L';( R,Z,m ) such that 


em(RZ = j fma (mV = Rẹ 9x VI) +97 (A41 


and 


E (RD =j (mY e Re q^ vh) x 9 | (A12 


This method of expanding the electromagnetic fields into modal 
components is named the Coupled Azimuthal Potential (CAP) method 
due to the fact that the scalar potentials T(R, Z, m) and T(R, Z, m) 
are coupled in these field equations and are not separable. 

Due to the symmetry in the azimuthal coordinate, 6, in radome 
design, an arbitrary meridional cross section of the radome 
construction may be represented in an (R, Z)-plane such that there 
exists a two-dimensional solution domain S, with a boundary curve 
forming the surface of revolution ôS. It has been shown by Morgan, 


Chang and Mei [Ref. 7] that in this solution domain the potentials, IT, 
and L5, satisfy the following formally self-adjoint pair of partial 


differential equations (PDE's), 


V. (f, Re, VI; )« «E; 2 mV.(£f,$xVD) 20 (A.1 3) 


and 


0 , (A.14) 


mV- (fox VR) + V.(£, Ry, VIz)+ +21 


with Dirichlet boundary conditions (BC) for IT, and I’, specified on òS. 
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In addition to this PDE formulation there is an easily managed 
variational principle (Ref. 7] established using either an 
Euler-Lagrange approach or a stationary-theorem approach. The 
variational principle is based upon the existence of a stationary 


functional of the potentials, I, and I5, and their gradients such that 


SF = || L(Rzn.r,vr.vr) dRaz, (A15) 
S 


where the Lagrangian operator, L, is unique within an arbitrary 
constant multiplier and arbitrary independent additive function. The 


Lagrangian operator has the quadratic form 


L=f, VI: (Re, VI + móxVI,) 


* f, VD: (Ru, V - mex VIj) 


l 2 2 
- z(e«m «uz ). (A16) 


With this analytical approach to the solution of the electromagnetic 
field's behavior through the axially symmetric radome structure, the 
determination of the field on a point-to-point basis is accomplished 


through a numerical finite-element algorithm. 


B. FINITE ELEMENT ALGORITHM 

The numerical solution of the electromagnetic field problem 
involving scattering by, and penetration of, the radome is 
accomplished by the development of a finite-element algorithm 
using the mesh geometry in Fig A.2. This configuration in the 


meridional (R, Z)-plane uses linear triangular elements which are 
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arranged by the EMRAD program to conform to the surfaces of the 
axisymmetric radome. The selection of surface conformation 
simplifies numerical analysis of the normal derivative discontinuities 
at the surface boundaries. The two-dimensional mesh configuration 
that is developed is named the semiannular conformal (SAC) mesh. 

Figure A.2 graphically represents the separate spatial regions 
which are of interest. The homogeneous inner spherical core of the 
radome occupies Region I. The inhomogeneous axisymmetric radome 
wall construction cross-section is in Region II. The unbounded 
exterior region exists exterior to the outer sphere in the 
homogeneous free-space of Region III. Note that the surfaces of the 
inner sphere, Region I, and the exterior Region III, bound the area of 
interest for numerical analysis in Region II. 

The electromagnetic fields within the homogeneous inner core 
of the radome (Region I) and the unbounded exterior (Region III) are 
defined in terms of the appropriate spherical harmonic modal 
expansions of the radial TE and TM field components. At the 
surfaces, rT = a and r = b, the modal elements in these fields are 
applied as sets of Dirichlet boundary conditions. Then, the numerical 
solution for the fields in the finite-element domain of Region II is 
generated for each boundary condition. Finally, the unknown modal 
expansion coefficients for the core and exterior scattered fields are 
obtained by enforcing continuity of the electromagnetic fields at the 


interfaces along r = r, andr=r,. This approach to numerical solution 


of the fields is the basis of the unimoment method. In the 


implementation of the variational principle using the finite-element 
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method, the azimuthal modal field components are selected as the 
fundamental unknowns. 

In Region II of Fig. A.2, the SAC mesh is composed of numerous 
interior nodes. Each of these interior nodes is surrounded by six 
triangular regions known as "elements". The azimuthal field 


components within the mesh are represented as basis function 


expansions, 
N 
¢-en(R,Z - Y e,(n) u,(RZ) (A.17) 
n=) 
and 
N 
ò- ha (R, = y ha (0) u (R2 . (A.1 8) 
n=l 


Each of the respective sets of e, (n) and h, (n) represent one of the N 
complex nodal values found in the azimuthal fields. The finite 


element basis functions u,(R, Z) are area coordinate linear pyramid 


functions which each assume a value of unity at their associated 
node, n, and zero at each of the neighboring nodes. There is an 
assumption of linear variation within each of the triangular elements 


in the surrounding group such that 


URA = aR Ebel Ca (A.1 9) 
The real constants (a, b, c) corresponding to each node, n, in element 
£ is found in terms of the nodal (R, Z)-coordinates of the element. 
The variational finite-element technique is initiated by 
Substituting the linear shape-function expansions in Equations (A.17) 


and (A.18) into the stationary Euler-Lagrange functional, SF, as 
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described in Equation (A.15). The stationary point of SF is then 
computed by imposing the null derivative condition with respect to 
each of the unknown coordinates representing the interior nodal 
fields, while constraining the "known" boundary nodal fields. The 


resulting linear system describing the nodal fields is 


N 
Y (es) J| [£, Re, V(Ru,):V(Ru;)— Re, u; JdRdZ 
i=l ST 


ch, GI (mf, V(Ru,)- 6x V(Ru;) dRdZ 20 (A.2.0) 
Sn 
and 


N 
S (hai) Íf [f£ Ru, V(Ru,). V(Ruj)- Riu, u; ] dRdZ 


i=] Sn 
~e,(i)| | [mf V(Ru,)- 6x V(Ruj)] dRdZ]-0 (A.2.1) 
Sn 
for each value of the internal nodes, n, where i spans all nodes, 
including boundary nodes. As these above equations are enforced 
for each internal node, they involve integrations over only the 


elements connected to this node, S,, and therefore relate the value of 
e,(n) and h,(n) only to adjacent nodal values, e,,(1i) and h, (1). Then, 


in the SAC mesh implementation, in the case of each internal node 
and the six adjacent nodes, the above equation will relate 14 
adjacent node azimuthal field values, some of which may be specified 


at the boundary. 
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The computation for each nodal value, as indicated by 
Equations (A.20) and (A.21), results in a sparse global system matrix 
which contains generally complex elements. At each node the mean 
values of the spatially variable e,(R, Z) and u,(R, Z) are used as the 
material constants in that triangular element. Applying the Dirichlet 
boundary conditions to the SAC mesh shown in Fig. A.2, Morgan and 
Mei [Ref. 5] found the numerical finite-element solution for the 
electromagnetic fields in Region II using a two-sweep block-by-block 
matrix inversion algorithm related to the Riccati transform for 
difference equations. At each constant colatitude, O = 60,, a vector of 


internal nodal azimuthal field values along the ray from r, to r, is 


defined by Bx. The resultant global system matrix defined by 
Equations (A.20) and (A.21) has a diagonal tri-block submatrix 
Structure, which can be represented locally by relationships between 


adjacent radial node vectors in the mesh, 


Ay Ba + E A +Q An = ak, (A.2 2) 
such that A,, B,, and C, are sparse submatrices. These form the 
tri-block diagonal structure of the global system matrix while ay, 
the local driving function, is generated from Equations (A.20) and 
(A.21) using the known Dirichlet BC's at r=a and r=b. 

The system solution in the two-dimensional meridional (R, Z) 
plane is performed in a step-wise manner beginning at the positive 
Z-axis (0, = 6, = O0) and proceeding in equal marching increments of 
A0. At each O,, the SAC mesh is constructed using a minimum 


distortion scheme to conform to the approximate piecewise smooth 
105 


surfaces. At each 8, increment, as the mesh is constructed in this 
sequential circumferential manner, the three block matrices are 
generated which relate adjacent radial node vectors. 

During the incremental marching of 0,, the block-by-block 
Riccati matrix inversion is implemented. The Riccati matrix, R,, and 
the associated vector, S,, are related to the vector, Bk, by the 
transformation 

Ba = Ri B, + $, . (A.2 3) 
Substituting this Riccati relation into Equation (A.22), solving for py 


and returning to the form of the Riccati transform yields the 


recurrence relations 


Ry = -( By + AyeR, ) 0G (A.24) 


and 


Sk = (Bi + AjeRi Ue Be -— Ajos, ). (A.2 5) 


Using these relations, the Riccati transform sweeps forward 
from initial conditions, using the equations of step-wise evolution in 
Equations (A.24) and (A.25) to sequentially generate and store the R, 
matrices and the S, vectors. This is called the "first sweep". The 


solution vectors B & are then obtained by using the transform in 
Equation (A.23) to perform a "backsweep" with calculations using the 
stored Riccati arrays. This sequentially localized technique for 


handling the Dirichlet boundary value problem blends easily with 
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the unimoment method which is used to solve the electromagnetic 


field in Regions I and III of Fig. A.2. 


C. UNIMOMENT METHOD 

The unimoment method is a computational procedure which 
provides for the numerical solution of the incident, scattered and 
spherical core electromagnetic field expansions when dealing with 
the constructed SAC finite-element mesh of Fig. A.2 The numerical 


solution is generated by representing the core (rsr,) and the exterior 
scattered (r2r,) fields in terms of radial TE and TM spherical 


harmonic expansions, each with unknown core and scattered modal 
field coefficients. In order to solve for these unknown coefficients, 
the core expansion of Region I is applied as a boundary condition 
along r = a, while the boundary condition at r = b is determined as 
the addition of the incident field and the scattered field expansion in 
Region III. 

Superposition is then applied to the boundary conditions, 
composed of the core field expansion and the exterior field expansion 
in order to find a complete solution. The complete field solution is 
found from the addition of the partial field solutions generated from 
two sets of applied BC's. The first set of boundary conditions applied 
is composed of the modal spherical harmonic core fields along r = a 
with zero fields along r = b. The second set of boundary conditions 
contains the spherical harmonic exterior scattered and incident fields 
at r = b with null fields at r = a. These multiple boundary conditions 


generate multiple solutions at the nodes along both boundaries, 
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r = T} adr = r>. The unknown core and scattered modal field 


coefficients are obtained by setting the weighted superposition of 
interior solutions equal to the original field expansions. This equality 


is enforced along the boundaries at r=r,andr=r,, utilizing least 


squares. 

The success of the unimoment method is contingent upon the 
capability to analytically represent the electromagnetic fields 
everywhere outside of the finite element mesh region. These 
representations, as developed for Regions I and III, are given as 
infinite series of spherical harmonics which must be truncated to 
facilitate numerical calculations. The plane-wave fields that are 
incident on the radome are decomposed in spherical coordinates into 


azimuthal modes as 


! "m jsin(m Q) 
E, (1,89) = 5 e s (59) (A.2 6) 
my cos(m 4) 
and 
| M | cos(m Q) 
n, H, (1,89) = X hi, ¿(1,0) ] (A.27) 
ss jsin (mg) 


where M is the mathematical truncation index of the series 
representation. The bracketed sinusoids represent the polarization 
convention used for TE and TM modes, respectively. 

During computations, the numerical solution of the unknown 
core and scattered modal field coefficients is run concurrently with 


the sequential Riccati transform finite-element solution in Region II. 
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On the forward sweep, boundary conditions for both the incident and 
the modal fields are generated only at the local node row being 
considered in that particular sweep step. On the backsweep the 
nodal values for each applied boundary condition are obtained along 


the r = r} and r = r, contours. The corresponding incident field and 


spherical harmonic field at each node are generated, and then the 
nodal residual fields are formed. During the backsweep, inner 
products integrations are accumulated and their contributions are 
added to a field moment matrix and the incident-field driving vector, 
which corresponds to each incident field. The field expansion 
coefficients for each incident field are obtained by multiplying the 
inverse of the field moment matrix with that field's driving vector. 
With the field expansion coefficients calculated and collected, the 
electromagnetic field within the radome core can be determined on a 


point-to-point basis as developed in Chapter II . 
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Appendix B 


RADSHP PROGRAM CODE 


PROGRAM RADSHP 


FOO GGG GGG GIG I i ake ak ak ak ak ak ak ak ak ak ak ad ld 


GENERATING THE MATERIAL FILE AND DATA STRUCTURE FILE 
FOR A LAYERED DIELECTRIC RADOME 
MATERIAL FILE PREPARED FOR EMRAD 
DATAFILE PREPARED FOR CURVE DIGITIZER 
USE EMCADIN TO GENERATE STRUCTURE FILE FOR EMRAD 
Written by R.J.Vince Nov-Dec 89 


ak k kc oc ke ok e dl ak al ae a e ae a a a a a a a a a al a al ok li ke KK KK KK KKK KK KK KK KK x 


VARIABLE DECLARATIONS 


ak e a ak ak ak ad a ak al a ok ae al a ae a ae a ke oco ae a Ik Ik k Ik ae a al a KK KKK Ik KK KKK KKK KKK KKK KKK KKK KK KKK 


O TCH O OOR GAO NHG QOO Oo 0 


REAL X(5,1001), Z(5,1001), PI, DELTHETA, DEGTORAD 

REAL THETA,C(5),R(5),H(5),DELZ(5),SHFT,FACTOR,F(5,7) 
INTEGER IBIG, I, L, SEPRATN, NSHAPE, NLMAX, RESOLV, RESPTS 
COMPLEX ER(5),UR(5) 

CHARACTER*10 SHAPE(6) 

CHARACTER*64 HDR1,HDR2 

CHARACTER*8 YNAME 

CHARACTER*12 MATFIL,DATFIL 

CHARACTER*1 YN 


$LARGE 


ak ak ak ak c kc kc ok ak ok ak oe OR kkk kk Rk kk kk 


INITIAL VALUES 


Add ll lalalala aaa alle ad ale al ad al ale ala ale al a a al ale ad af ol al af al ol al a a ol ol 


@ aC) oO) 
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PI=3.1415927 
DEGTORAD = PI/180. 
SEPRATN = 999990 
RESOLV=1000 
RESPTS=RESOLV+1 
NLMAX = 5 


C ITITIALIZE SHAPE FUNCTION FACTOR ARRAY 


E 
DO 67 I=1,NLMAX 
DO 65 J=1,6 
F(1,)20.0 
65 CONTINUE 
67 CONTINUE 
C 
C 
(Rad OIGO ak 2e ake ake ak 3k ak ak ak ak ak ak ak akak ak ak ake ak akak ak ak 
C INPUT NECESSARY DATA 
Co ess aga oda a ak ak ak ak ak ak ak a ak ak ak ake akak ak ak ak ak akak sk sk 3k 3k ske ake ake ak ake ake ake ak ak ak ak ake ak 3k Sk ak ak ak akak ak ak ae ak akak ak k 
C 
C PROMPT FOR OUTPUT FILENAMES 
C 
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pd 


QN mn Bb W N 


WRITE(*,110) 
FORMAT(/9X,'WELCOME TO RADSHP", 
//9X, THIS PROGRAM GENERATES RADOME MATERIAL AND SHAPE DATA’, 
//(7X FILES FOR USE WITH EMRAD, BASED ON USER INPUT.', 
/I//9X, TWO FILES ARE REQUIRED TO HOLD PROGRAM OUPUT DATA.., 
JINX, PLEASE ENTER THE NAME OF THE FILE TO BE USED TO ', 
[TX HOLD THE PROGRAM MATERIAL PARAMETER DATA. THE EXTENSION ', 
/X,.MAT WILL BE ADDED AUTOMATICALLY, I. E. YOURNAME.MAT'/) 
READ(*,100) YNAME 
NCHAR - INDEX(YNAME; ) -1 
MATFIL = YNAME(1:NCHAR) // '.MAT 
OPEN(1,FILE-MATFIL,STATUS-'UNKNOWN', ACCESS-'SEQUENTIAL', 
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111 


C 


FORM-'FORMATTED) 
WRITE(*,111) 

FORMAT(///9X,'PLEASE ENTER THE NAME OF THE FILE TO BE USED ', 
/1X,'TO HOLD THE PROGRAM OUTPUT DATA. THE EXTENSION .DAT WILL’, 
/TX,'WILL BE ADDED AUTOMATICALLY, I. E. YOURNAME.DAT /7X) 
READ(*,100) YNAME 

NCHAR = INDEX(YNAME;  -1 

DATFIL » YNAME(I:NCHAR) // .DAT 
OPEN(2,FILE=DATFIL,STATUS='UNKNOWN' ACCESS='SEQUENTIAL,, 
FORM="FORMATTED’) 


C PROMPT FOR # OF LAYERS WITHIN THE RADOME 


C 
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115 


C 


WRITE(*,114)NLMAX 
FORMAT(//X,ENTER THE NUMBER OF LAYERS (.LE. 'I2,) //7X) 
READ(*,*) NL 


IF ( (NL .LT. 1) .OR. (NL .GT. NLMAX) ) THEN 
WRITE(*,115)NLMAX+1 

FORMAT(//7X,'NUMBER OF LAYERS IS OUT OF RANGE, (0<#<',12,)'//7X) 

GOTO 5 

ENDIF 


C LOOP TO INPUT LAYER INFORMATION 


C 


tu d» ù N 


WRITE(*,7) 
FORMAT(/8X,Each layer of the radome will have material consts,’, 
/10X,'an outer height of the radome layer,H,’, 
/10X,'a outer base radius,R, and’, 
/10X,'a shape function, & shape factor parameters (if appl.),', 
//10X,'The first "layer" will be the inner core of the radome, ”, 
/10X,'normally filled with air.’) 
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39 


DO 12 L=1,NL 

CONTINUE 

WRITE(* 9) L 

FORMAT(/,1 Layer No: ',13) 


IF (L.EQ.1) THEN 

WRITE(*,*)' . This is the inner core of the radome' 
ELSEIF (L.EQ.NL) THEN 

WRITE(*,*) Thisis the outer layer of the radome' 
ENDIF 


C Prompt for Height and Radius of this layer 


C 


43 


WRITE(*,*) 
WRITE(*,43) 

FORMAT(/8X, Please enter the radome layer height, H,' 
/8X, and the outer base radius,R, for this layer.) 
READ(*,*) H(L),R(L) 


C Prompt for shaping function, to determine equations to be used 


€ 
(o 


10 


to generate this layer of the radome 


oY KN tA A WY N 


WRITE(*,*) 

WRITE(*,10) 

FORMAT(/8X, The following shaping functions are available, ' 
/8X, | 1. General Ogive' 


/(8X, 2. Tangent Ogive' 
/8X, 3. von Karman’ 
/8X, | 4. Power Series 
/8X,' | 5. Parabola, curved' 


/8X,' 6. Parabola, capped by pointed apex' 
//8X,  Indicate desired shaping function by entering a number' 
/MOX, 132, ... or 6) 

S 


READ(*,*) NSHAPE 


IF (NSHAPE.EQ.1) THEN 
Ogive 
SHAPE(L) = 'Gen. Ogive' 
CALL OGIVE(L,X,Z,H,R,F,RESOLV,RESPTS NLMAX) 


ELSEIF (NSHAPE.EQ.2) THEN 
Tangent Ogive 
SHAPE(L) = "TAN Ogive ' 
CALL TANOGIVE(L,X,Z,H,R,F,RESOLV ,RESPTS,NLMAX) 


ELSEIF (NSHAPE.EQ.3) THEN 
von Karman 
SHAPE(L) = 'von Karman’ 
CALL VONKRMN(L,X,Z,H,R,RESOLV,RESPTS, NLMAX) 


ELSEIF (NSHAPE.EQ.4) THEN 
Power Series 
SHAPE(L) = Pwr Series’ 
CALL POWER(L,X,Z, H,R,F,RESOLV,RESPTS,NLMA X) 


ELSEIF (NSHAPE.EQ.5) THEN 
Parabola (curved) 
SHAPE(L) = 'Para Curve' 
CALL PARACRV(L,X,Z,H,R,F,RESOLV,RESPTS NLMAX) 


ELSEIF (NSHAPE.EQ.6) THEN 
Parabola with apex 
SHAPE(L) = ‘Para APEX’ 
CALL APEX(L,X,Z,H,R,F,RESOLV,RESPTS NLMAX) 


ELSE 
Selection out of limits 
WRITE (*,*) ‘Selection is not available, try again.’ 
WRITE (*,*) 
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GOTO 39 
ENDIF 


C Provide listing of inner layer parameters 


C 


11 


12 


C 


WRITE(*,554) 
DO 11 Iz1,L 

WRITE (*,555) LSHAPE(D,H(D,R(I)(F(1J)J21,7) 
CONTINUE 


WRITE(*,*) 


CONTINUE 


C SHIFT RADOME SO THAT ORIGIN IS NEAR THE CENTER OF INNER LAYER 
C NOTE: FACTOR is adjusted for optimal placement of origin for 


C 
C 


mesh construction by EMRAD 


FACTOR=0.4 


C Max inner height = Z(NL,RESPTS) 


C 


17 
88 


OEEO 


SHFT=Z(1,RESPTS)*(FACTOR) 


DO 88 L=1,NL 
DO 77 I=1,RESPTS 
Z(L,)-Z(L.I)-SHFT 
CONTINUE 
CONTINUE 


3k kc okc ok okc oc oc okc oc okc oc ok ok KKK HHH KK HK KKK HK HK KKH KKK KKK KKK KKH KH KKK HK HK a ad ok oe ok ok oko oe o 


OUTPUT 


3c kc kc e a ok okc ok e ae ad okc oc okc okc ok okc ok ok al ae al add ae al e al ad ad ad ad ad oc ok ad a ad ad de a ad ad ad ad oc ok a ad oc ok ok ok oe oe oe oe ok oe o o 9 
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DO 330 L=1,NL 
C OUTPUT LAYER X,Z POINTS FROM TIP TO X-Z PLANE 
DO 220 I=RESPTS,1,-1 
WRITE(2,*) X(L,D,UL,I) 
220 CONTINUE 


C CLOSE LAYER BELOW SHIFTED X-Y PLANE 


C  ATZ=[-1*[SUM OF (LAYER THICKNESS)] - SHFT ] 
C 
E 
C DETERMINE LAYER THICKNESS, INNER LAYER HAS NO SHIFT FOR THICKNESS 
C OTHER LAYERS HAVE BOTTOM THICKNESS DETERMINED BY RADIAL THICKNESS 
C NOTE - THE SHIFT PERFORMED FOR THICKNESS IS ADDITIVE FOR LAYERS 
E 
IF (L.EQ.1) THEN 
THICK =0 
ELSE 
THICK = THICK + R(L)-R(L-1) 
ENDIF 
$ 


BOTTM - (-1.0 * THICK) - SHFT 

WRITE(2,*) R(L), BOTTM 

WRITE(2,*) 0.0, BOTTM 

WRITE(2,*) SEPRATN,SEPRATN 
330 CONTINUE 


C SIGNAL TERMINATION OF DATA SET 
WRITE(2,100) "END, END' 


C Put summary of data to screen 
C 
WRITE(*,554) 
DO 440 I=1,NL 
WRITE (*,555) I,SHAPEOD,HO),RO),(F(1,) J=1,7) 
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440 CONTINUE 


C 
C PROMPT FOR HEADER INFORMATION 
E 
WRITE(*,112) 
112 FORMAT(//7X, THE HEADERS ALLOW THE USER TO IDENTIFY THIS'//7X, 
1'SET OF DATA FROM ALL OTHER SETS.., 
2//7X PLEASE ENTER HEADER +1 (64 CHARACTERS MAX) //7X) 
READ(*,100) HDR1 
WRITE(*,113) 
113 FORMAT(/7X,PLEASE ENTER HEADER +2 (64 CHARACTERS MAX)'//7X) 
READ(*,100) HDR2 


O 


WRITE HEADER INFORMATION INTO MATERIAL OUTPUT DATA FILE 


WRITE(1,105) HDR1 
WRITE(1,105) HDR2 


O 


WRITE LAYER INFORMATION INTO MATERIAL OUTPUT DATA FILE 


WRITE(1,106) NL 


O 


Prompt for material Er for each layer of construction 


DO 812 L=1,NL 
WRITE(*,9) L 


IF (L.EQ.1) THEN 
WRITE(*,*)' This is the inner core of the radome' 
ELSEIF (L.EQ.NL) THEN 
WRITE(*,*)' This is the outer layer of the radome 
ENDIF 
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C Prompt for material constants 


WRITE(*,*) ‘Enter Er for this Layer: ' 
READ (*,*) ERR 


C Ensure that Imag Er is < 0 


ERI = -0.000001 


ER(L)=CMPLX(ERR.ERI) 

UR(L)=CMPLX(1.0,-1.E-6) 

WRITE(1,102) ER(L),UR(L) 
812 CONTINUE 


C PLACE COPY OF CONSTRUCTION AT BOTTOM OF MATERIAL FILE 
E 
WRITE(1,554) 
DO 660 I=1,NL 
WRITE (1,555) I,SHAPE(1),H(),R(),(FCJ) J=1,7) 
660 CONTINUE 


C --- FINAL OUTPUT TO SCREEN 


WRITE(*,999) MATFIL,DATFIL,DATFIL 

999 FORMATY/8X, 
1' Computations and output are completed,the datafile ,',A, 
2/,8X, holds the material specifications for input to EMCAD.' 
3//8X, The datafile, ',A,', holds the radome shape data for 
4/,8X, viewing in the CAD package, CURVE DIGITIZER. ' 
5//,8X, The structure file required for input to EMCAD can be' 
6/,8X, generated by the executing the program EMCADIN, 
7/,8X, using, ',A,', as input to EMCADIN, ' 
8//,8X,' The output file from EMCADIN will be in proper format 
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9/,8X,' for input to EMCAD for the structure file.") 
WRITE(*,*) 
C 
Ç alakaa aaao aoaaa oook 
100 FORMAT(A) 
101 FORMAT(15) 
102 FORMAT(4(E14.6)) 
103 FORMAT(IS,S(E14.6)) 
104 FORMAT(2(15)) 
105 FORMAT( ',A) 
106 FORMAT(IS) 
554 FORMAT (Q//IIIIIIMIIIII 
1 //1X,LAYER',2X,' SHAPE ',2X, HEIGHT',2X,RADIUS', 
EZX EBI 2X, F2 2X, F3'2X, FA'2X, FS ', 
22X, F6'2X, F7 '/1X, ----- 2X, ---------- 2X, ------ ña 


A DD X >) 


555 FORMAT (1X,15,2X,A10,2X,F6.2,2X,F6.2,7(F6.2)) 


9999 CONTINUE 

C 

C END PROGRAM RADOME 

E 
CLOSE(1) 
CLOSE(2) 
STOP 
END 

RARA 

C OGIVE SHAPING FUNCTION SUBROUTINE 
SUBROUTINE OGIVE(L,X,Z,H.R,F,RESOLV,RESPTS,NLMAX) 
INTEGER IBIG, I, L, SEPRATN, NSHAPE, NLMAX, RESPTS, RESOLV 
REAL X(NLMAX RESPTS), Z(NLMAX,RESPTS), PI, DELTHETA, DEGTORAD 
REAL THETA,C(5),R(5),H(5),DELZ(5),SHFT,FACTOR 
REAL ALPHA1(5), ALPHA(5), B, M, RLIGAMMA, F(NLMAX,7) 


COMPLEX ER,UR 
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CHARACTER*1 YN 


C SET CONSTANTS 


PI=3.1415927 
DEGTORAD = PI/180. 


P 


DEFINE START AND END POINT OF EACH RADOME LAYER 


Z(L,1)=0.0 
X(L,1)=R(L) 
Z(L,RESPTS)=H(L) 
X(L,RESPTS)=0.0 


s DEFINE RADOME LAYER BY OGIVE EQUATIONS .... 


DELZ(L)-H(L)/RESOLV 


C Prompt for Shaping factors for this ogive layer 
C 
WRITE(*,*) 
WRITE(*,43) 
43 FORMAT(8X, The ogive shaping function requires three shaping’, 
1/8X factors for determining the radome construction.’ 
1//8X,'The first factor, F1, determine the arc of curvature’, 
1/8X,'of the radome surface. Please enter F1:") 
READ(*,*)F(L,1) 


WRITE(*,*) 
WRITE(*,45) 
45 FORMAT(8X,' Two other shaping factors are used to determine the’, 
1/8X , relative position of the center of the ogive shaping’, 
1/8X,'arc, relative to the radome apex.', 
1//8X,' The first of these factors, F2, determines the perp. ', 


1/8X, distance from the z-axis to the center of the arc.', 
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1//8X, The second positioning factor, F3, determines the’, 

1/8X,the distance from the center of the arc to the apex’, 

1/8X,'of the radome, along the z-axis.', 

1//8X,'Please enter the position factors for the center of the’, 

1/8X, curvature arc, F2 and F3. 
READ(*,*)F(L,2),F(L,3) 


DO 33 I=2,RESPTS-1 


Z(L,D=Z(L,I-1)+DELZ(L) 


C EQUATION OF RADOME IS x-SQRT( F1**2-(F2-(H-z))**2 ) - F3 


X(L,I) = SQRT(F(L,1)**2 -(F(L,2) - (H(L)-Z(L,D) )**2 )-F(L,3) 


33 CONTINUE 


C ...... END OF RADOME LAYER .... 
RETURN 
END 

jaooo END OF OGIVE 

a 

C TANGENT OGIVE SHAPING FUNCTION SUBROUTINE 
SUBROUTINE TANOGIVE(L,X.Z,H,R,F,RESOLV,RESPTS,NLMAX) 
INTEGER IBIG, I, L, SEPRATN, NSHAPE, NLMAX, RESPTS, RESOLV 
REAL X(NLMAX,RESPTS), Z(NLMAX,RESPTS), PI, DELTHETA, DEGTORAD 
REAL THETA,C(5),R(5),H(5),DELZ(5),SHFT,FACTOR 
REAL ALPHA 1(5), ALPHA(5), B, M, RL,GAMMA, F(NLMAX,7) 
COMPLEX ER,UR 
CHARACTER*1 YN 


C SETCONSTANTS 


PI=3.1415927 


DEGTORAD - PI/180. 
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C 


C DEFINE START AND END POINT OF EACH RADOME LAYER 


C 


Z(L,1)=0.0 
X(L,1)=R(L) 
Z(L,RESPTS)=H(L) 
X(L,RESPTS)=0.0 


.... DEFINE RADOME LAYER BY TANGENT OGIVE EQUATIONS .... 


DELZ(L)=H(L)/RESOLV 


Determine Shaping factors for this ogive layer 


F(L,4) = R(L)/2.0 + 2.0*( H(L)**2.0 /R(L) 
F(L,5) = F(L,4) - R(L)/4.0 


WRITE(*,*) 
WRITE(*,43)F(L,4),F(L,5) 


43 FORMAT(8X, The tangent ogive shaping function determines two’, 
1/8X, ‘factors for the radome construction. These are calculated’ 
1/8X, based on the height and radius of the radome layer.', 
1/8X,'These factors are F4 = ',F6.2,'and F5 = ',F6.2) 


DO 33 I=2,RESPTS-1 


Z(L,D=Z(L,I-1)+DELZ(L) 


EQUATION OF RADOME IS x=SQRT( F4**2 - z**2 ) - FS 


X(L,I) 2 4.0* (SQORT( F(L,4)**2- Z(LjD**2 ) - FL,5)) 


33 CONTINUE 


.... END OF RADOME LAYER .... 
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RETURN 
END 
a SD) OF TAN OGIVE 
an 
C von Karman SHAPING FUNCTION SUBROUTINE 
SUBROUTINE VONKRMN(L,X.Z,H,R,RESOLV,RESPTS NLMAX) 
INTEGER IBIG, I, L, SEPRATN, NSHAF:-, NLMAX, RESPTS, RESOLV 
REAL X(NLMAX,RESPTS), Z(NLMAX RESPTS), PI, DELTHETA, DEGTORAD 
REAL THETA,C(5),R(5),H(5),DELZ(5), SHFT,FACTOR,SQP,ZF 
REAL ALPHA1(5), ALPHA(5), B, M, RL,GAMMA 
COMPLEX ER,UR 
CHARACTER*1 YN 


C SET CONSTANTS 
PI=3.1415927 
SQP=SQRT(PI) 
DEGTORAD = PI/180. 
C DEFINE START AND END POINT OF EACH RADOME LAYER 
Z(L,1)=0.0 
X(L,1)=R(L) 
Z(L,RESPTS)=H(L) 
X(L,RESPTS)=0.0 
C ...... DEFINE RADOME LAYER BY VON KARMAN EQUATIONS .... 
DELZ(L)=H(L)/RESOLV 
DO 33 I=2,RESPTS-1 
Z(L,D=Z(L,I-1)+DELZ(L) 
C EQUATION OF RADOME IS 


E x=(R/(SQRT(PI) ) * SQRT( FZ - SIN(2*FZ)/2 ) 
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C for FZ = ACOS( 1 - 2 * (H-z)/H) 
FZ =ACOS(1.0-2.0 *( H(L)-Z(L,I) )/H()) 
X(L,D) = R(L)/SQP * SQRT( FZ - 0.5*SIN(2*FZ) ) 


33 CONTINUE 


C ...... END OF RADOME LAYER .... 
RETURN 
END 

ak ak ak ak ak akk akk ae akk ak ak akk akk ak ak ak ak k k k k E NIT) OF VON KARMAN 

EA ARO 

C POWER SERIES SHAPING FUNCTION SUBROUTINE 
SUBROUTINE POWER(L,X,Z,H.R,F,RESOLV,RESPTS,NLMAX) 
INTEGER IBIG, I, L, SEPRATN, NSHAPE, NLMAX, RESPTS, RESOLV 
REAL X(NLMAX,RESPTS), Z(NLMAX RESPTS), PI, DELTHETA, DEGTORAD 
REAL THETA,C(5),R(5),H(5), DELZ(5),SHFT,FACTOR 
REAL ALPHA 1(5), ALPHA(5), B, M, RL,GAMMA, F(NLMAX,7) 
COMPLEX ER,UR 
CHARACTER*1 YN 


C SET CONSTANTS 


PI=3.1415927 
DEGTORAD = PI/180. 


DEFINE START AND END POINT OF EACH RADOME LAYER 


O 


Z(L,1)=0.0 
X(L,1)=R(L) 
Z(L,RESPTS)=H(L) 
X(L,RESPTS)=0.0 


Coss DEFINE RADOME LAYER BY OGIVE EQUATIONS .... 
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DELZ(L)=H(L)/RESOLV 


C Prompt for Shaping factor for this power series layer 
C 

WRITE(*,*) 

WRITE(* 43) 

43 FORMAT(8X, The power series shaping function requires a', 
1/8X,'shaping factor, F6, to be used in the exponential in’, 
1/8X,'order to determine the radome construction." 
1//5X, When F6 = 1.0 the radome shape is conical’, 
1/8X, For F6 = 0.75 the radome shape approximates the ', 
1/8X,'Newtonian contour (minimizes drag). Please enter F6:’) 

READ(*,*) F(L,6) 


DO 33 I=2,RESPTS-1 


2(L,D=Z(L,I-1)+DELZ(L) 


C EQUATION OF RADOME IS x-R * ( (H-z)/H )**F6 


X(L,I) 2 R(L) * ((H(L)-Z(L,D) / H(L) )**F(L,6) 


33 CONTINUE 


C... END OF RADOME LAYER .... 
RETURN 
END 
exóoooooo|teoeeeeee END OF POWER SERIES 
RENE eoa inane oaa oce 
C PARABOLIC SHAPING FUNCTION SUBROUTINE 
SUBROUTINE PARACRV(L,X,Z,H,R,F,RESOLV,RESPTS NLMAX) 
INTEGER IBIG, I, L, SEPRATN, NSHAPE, NLMAX, RESPTS, RESOLV 
REAL X(NLMAX RESPTS), Z(NLMAX RESPTS), PI, DELTHETA, DEGTORAD 
REAL THETA,C(5),R(5),H(5),DELZ(5),SHFT,FACTOR 
REAL ALPHA1(5), ALPHA(5), B, M, RL,GAMMA, F(NLMAX;,7) 
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COMPLEX ER,UR 
CHARACTER*1 YN 


SET CONSTANTS 


PI=3.1415927 
DEGTORAD = PI/180. 


DEFINE START AND END POINT OF EACH RADOME LAYER 


Z(L,1)=0.0 
X(L,1)=R(L) 
Z(L,RESPTS)=H(L) 
X(L,RESPTS)=0.0 


M DEFINE RADOME LAYER BY PARABOLIC EQUATION .... 


Prompt for Shaping factor for this parabolic layer 


WRITE(*,*) 
WRITE(*,43) 

43 FORMAT(8X,' The parabolic shaping function requires a', 
]/8X,'shaping factor, F7, to determine the curvature of, 
1/8X,'the parabolic arc for the radome construction.’ 
1//8X Please enter F7:') 

READ(*,*)F(L,7) 


iem FOR ROUNDED PARABOLIC RADOME LAYER .... 


... RADOME IS SCALED TO A HEIGHT = SCALE FUNCTION,C 


DELZ(L)=F(L,7)/RESOLV 


DO 33 I22,RESPTS-1 
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Z(L,N=Z(L,]-1)+DELZ(L) 


C RECALL EQUATION OF RADOME IS X^2=F7-Z, X=SQRT(F7-Z) 
X(L,D=SQRT(F(L,7)-Z(L,D ) 
33 CONTINUE 
C 
C SCALE POINTS TO RADOME DIMENSIONS 
C 


DO 44 I=2,RESPTS-1 
X(L,D)=X(L, I)*R(L/SORT(F(L,7) ) 
Z(L.D-Z(,D* HOLVEQ,7) 
44 CONTINUE 


gs END OF ROUNDED RADOME LAYER .... 


ao: END OF RADOME LAYER .... 
RETURN 
END 

aoaaa END OF PARA CURVE 

k kckkolokekekoolokekekokookokolokololokokolololololololokolololek 

C POINTED PARABOLIC SHAPING FUNCTION SUBROUTINE 
SUBROUTINE APEX(L,X,Z,H,R,F,RESOLV,RESPTS NLMAX) 
INTEGER IBIG, I, L, SEPRATN, NSHAPE, NLMAX, RESPTS, RESOLV 
REAL X(NLMAX,RESPTS), ZINLMAX,RESPTS), PI, DELTHETA, DEGTORAD 
REAL THETA,C(5),R(S)H(S),DELZ(S),SHFT,FACTOR 
REAL ALPHA 1(5), ALPHA(5), B, M, RL,GAMMA, F(NLMAX, 7) 
COMPLEX ER,UR 
CHARACTER*1 YN 


C SET CONSTANTS 


C 
M27 


PI=3.1415927 
DEGTORAD = P1/180. 
E 
C DEFINE START AND END POINT OF EACH RADOME LAYER 
C 
Z(L,1)=0.0 
X(L,1)=R(L) 
Z(L,RESPTS)=H(L) 
X(L,RESPTS)=0.0 


Cs DEFINE RADOME LAYER BY POINTED PARABOLIC EQUATION .... 


C Prompt for Shaping factor for this parabolic layer 


WRITE(*,*) 
WRITE(*,43) 
43 FORMAT(8X,' The APEX parabolic shaping function requires a’, 
1/8X 'shaping factor, F7, to determine the curvature of, 
1/8X,'the parabolic arc for the radome construction and the’, 
1/8X, percent of the radome that is the pointed apex.’, 
1//8X,Please enter F7:') 
READ(*,*)F(L,7) 


C ...RADOME IS SCALED TO A HEIGHT = MODIFIED SCALE FUNCTION, F7+1 
DELZ(L)=( F(L,7)+1.0 /RESOLV 


DO 55 I=2,RESPTS-1 


Z(L,D=Z(L.I-1)+DELZ(L) 


IF( Z(L,D.LE.( F(L,7) -1 ) ) THEN 


C RECALL EQUATION OF CONE IS X^2=F7-Z, X=SQRT(F7-Z) 
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X(L,D-SQRT( F(L,7)-Z(L,I) ) 
ELSE 
C STRAIGHT LINE FORMULA FOR POINTED END 
X(L,D=( (F(L,7) +1 3-Z(L,D )/2.0 
ENDIF 


55 CONTINUE 


C 
C SCALE POINTS TO RADOME DIMENSIONS 
C 
DO 66 IZ2,RESPTS-1 
X(L,D-X(L,D*R(L)/SQRT( F(L;7) ) 
ZL,D=Z(L,D*H(L)A F(L,7)+1.0 ) 
66 CONTINUE 


Ce... END OF POINTED RADOME LAYER .... 


Cua END OF RADOME LAYER .... 


END 
JOGO KK END) OF PARA APEX 
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Appendix C 


CORFLD PROGRAM CODE 
PROGRAM CORFLD 


Program Created by LT R J. Vince (Aug-Dec 89) 
Modifications by M.A. Morgan (Nov-Dec 89) 


Computing core E-field at each point along a 50x50 radome grid 
in the radome's field coordinate system. 


Ore, Ce) 


INTEGER BELL,FF,NPTS,L,J,K,K1,NA,NINC 
INTEGER MSTART,MSTOP,N1,NM1,M,MF,MD,M1,M2,MM 
INTEGER I,ITE,ITM,KN,CA,CB,NJN,NLP,SELECT,SELCT1 
REAL ALPHA(13),DALPHA(13),ANGINC(13) 
REAL THETA,PHI,PI,DTR,RMIN,RCORE 
REAL XP,YP,ZP,R,COST,SINT,COSP,SINP,S,C,CMP,SMP,SA,SAMX 
REAL R1,R2,RSQ,EM,TH,FI,RPTS,L1,J1,PLOT(50,50),A(3,3) 
REAL P(40),DP(40),TP(6),PP(6) 
REAL EOTM.EOTE,EI(8),EDTERR(50,50,2), EDTER1,EDTER2 
COMPLEX EXP1,E(8),EDOT(S0,50,2) EDOT1,EDOT2 
COMPLEX EP1,EP2,ET1,ET2,ER1,ER2 JAY,RK,SRE,SRU,KR 
COMPLEX ERAD/(50,50,2),ETH(50,50,2),EPHI(50,50,2) 
COMPLEX ER(16),UR(16),COEF(70,80),CJ(40),DCJ(40), CPLOT 
CHARACTER*1 YN 
CHARACTER*8 PRTDAT 
CHARACTER*12 PRTDT1,PRTDT2,PRTDT3 
CHARACTER*64 GRFLAB,HDR1,HDR2,TITLE1,TITLE2,DUM 
$LARGE 
C 
C Setting constants 
BELL=CHAR(7) 
FF = CHAR(12) 
PI = 3.14159265359 
DTR= PI/180. 
JAY= (0.0,1.0) 


C Set default angle of rotation of planar antenna array away from z axis 


C Set default aspect ratio for plotung 
ASPECT = 1.35 


C Set number of points of resolution across array 
NPTS=50 
RPTS=( FLOAT(NPTS) - 1. )/2. 


C** Initialize arrays 
C Initialize the E-field matrices 
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C 
DO 1 L=1,NPTS 
DO 1 J=1,NPTS 
DO 1 K=1,2 
ERAD(J,L,K)=(0.,0.) 
ETH(,L,K)=(0.,0.) 
EPHI(J,L,K)=(0.,0.) 
1 CONTINUE 


C Prompt for input data filename without extension 
C 
WRITE(*,14) 
READ(*,101) PRTDAT 
NCHAR-INDEX(PRTDAT, ? -1 
PRTDT1 = PRTDAT(1:NCHAR) // '.INT 
OPEN(4 FILE=PRTDT1) 


TITLE1=The datafile used to generate fields is ' 
TITLE1=TITLE1(1:40) // PRTDT1 


PRTDT2= PRTDAT(1:NCHAR) //'.DAT' 
OPEN(15,FILE=PRTDT2) 


C Write header from structure data file 

e 
READ(14,101) HDR1 
READ(14,101) HDR2 
READ(14,101) GRFLAB 
WRITE(15,101) HDR1 
WRITE(15,101) HDR2 
WRITE(15,101) GRFLAB 
WRITE(*,101) HDR1 
WRITE(*,101) HDR2 
WRITE(*,101) GRFLAB 


C Read interior core Er,Ur 

C 
READ(14,507) ER(1),DUM 
READ(14,507) UR(1),DUM 
SRE=CSQRT(ER(1)) 
SRUZCSQRT(UR(1)) 
KR=SRE*SRU 
WRITE(*,*) 
VWaARTTE(*;*)'SRE= ';SRE 
WRITE(*,*)'SRUz ;SRU 
WRITE(*,*yKR =' KR 
WRITE(15,*)' SREZ,SRE,'; SRU=',SRU,'; KR = ',KR 


C Read minimum diameter of missile radome core expansion 
C Assuming in Ko*r Wavenumber Units 
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C 
READ(14,503) RMIN,DUM 
RCORE-RMIN 


C Read number of incident field angles (.LE. 5) 
C 

READ(14,500) NA,DUM 

NINC-2*NA 
C Read incident angles 

IPZ=0 

IMZ=0 

DO 11 I=1,NA 

READ(14,511) DANGLE 

C Checking for +/- Z-Axis Incident Fields 


IF(DANGLE.EQ.0.0) IPZ=I 
IF(DANGLE.EQ.180.0) IMZ=I 


DALPHA(I)= DANGLE 


EMRAD uses Poynting vector angle with Z-axis reference 
ANGle of incidence = 180 - DALPHA from EMRAD 


C050 


ANGINC(I)- 180- DANGLE 


WRITE(*,*) ' Incoming Inc. Angle #',] 

WRITE(*,*) ‘'EMRAD Poynting angle of incidence = ',DALPHA(I) 
WRITE(*,*) ' Radome incident angle = ' ANGINC(I) 
WRITE(*,*) 

WRITE(15,*) ' Incoming Inc. Angle 2I, 2 ,DALPHA(I) 


C Set-up array of angles in radian measure 
E 
ALPHA(D)=DTR*DALPHA(I) 
11 CONTINUE 
NAI=1 
IF (NA.EQ.1) GO TO 13 


C Determine angle of interest 
C 
12 CONTINUE 

WRITE(*,*) ' Select Angle of Interest By Number ' 
WRITE(*,*) 
READ(*,*) NAI 
IF ((NALGE.1).OR.(NALLE.NA)) GO TO 13 
WRITE(*,*) ' That Number is Not Available ... Try Again’ 
WRITE(*,*) 
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GO TO 12 


13 CONTINUE 


C 
E 


C 
C 


Checking if Selected Field is +/- Z-Axis Incident 
IM1=0 
IF((NAT.EQ.IPZ).OR.(NAI.EQ.IMZ)) IM1=1 


WRITE(*,*) 

WRITE(*,*) ' The angle of consideration is angle +" NAI 
1 , =',ANGINC(NAI),’ degrees’ 

WRITE(*,*) 


WRITE(15,*) ' The angle of consideration is angle #',NAI 
1 , =',ANGINC(NAI),’ degrees’ 
WRITE(15,*) 


Transformation of coordinates to gimble mounted planar array 


WRITE(*,*)' The planar array is set by default to lie in the’ 
WRITE(*,*) x-y plane, with boresight in the z-direction.’ 
WRITE(*,*) 
WRITE(*,*)' Do you want to change the array boresight ?' 
WRITE(*,*) ("Y" or "N")' 
IFLG=0 
DTH=0.0 
READ (*,101) YN 
IF(YN.EQ.'Y'".OR.YN.EQ.'y) THEN 
WRITE(*,*) 
WRITE(*,*) 'Enter Boresight Theta Angle (in degrees)’ 
READ(*,*) DTH 
TH = DTH * DTR 
WRITE(*,*) 
WRITE(*,*) ‘Enter Boresight Phi Angle (in degrees) 
READ(*,*) DFI 
FI 2 DFI * DTR 


IFLG=1 
Loading Transformation Matrix from Array to Core Coordinates 


CT=COS(TH) 

ST=SIN(TH) 

CP=COS(FI) 

SP=SIN(FI) 
A(1,1)=CT*CP**2.0+SP**2.0 
A(1,2)=CP*SP*(CT-1.0) 
A(1,3)=ST*CP 
A(2,1)=CP*SP*(CT-1.0) 
A(2,2)=CT*SP**2.0+CP**2.0 
A(2,3)=ST*SP 
A(3,1)=-ST*CP 
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A(3,2)=-ST*SP 
A(3,3)=CT 


ENDIF 


C Get modal information 


e 


READ(14,505) MSTART,MSTOP,DUM 
WRITE(*,*) MSTARTZ,MSTART, MSTOPZ',MSTOP 
WRITE(15,*)' MSTARTZ,MSTART, MSTOPZ',MSTOP 


C Loop modals to read modal coefficients 


C 


WRITE(*,*) 'Reading Coefficients and Constructing 
WRITE(* *) Fields for Each Azimuthal "m" Mode' 
WRITE(*,*) 


IF((MSTART.EQ.0).AND.(IM1.EQ.1)) THEN 
WRITE(*,*) ‘Terminating M-Loop at m=1 for +/- Z-Axis Incidence’ 
WRITE(*,*) 
MSTOP=1 
ENDIF 


DO 777 MD=MSTART,MSTOP 


READ(14,500) M,DUM 

WRITE(15,*)' M=',M 

IF(M.EQ.MD) GO TO 220 

WRITE(*,*) 'Error in Reading "m" - Check Data File’ 
STOP 


220 CONTINUE 


WRITE(*,*) 
WRITE(*,*) 'Azimuthal Mode "m" =',M 
WRITE(* *) 
IF(M.EQ.0) THEN 
EM= 1.0 
MF = 1 
ELSE 
EM = 2.0 
MF = 0 
ENDIF 


WRITE(*,*) ‘Reading Coefficients’ 

READ(14,500) NM1,DUM 

DO 665 ID=1,NA 

READ(14,500) IDUM 

IF(ID.EQ.D GO TO 221 

WRITE(*,*) 'Inc Angle Index Read Error - Check Data File' 
STOP 


221 CONTINUE 


WRITE(15,*) Iz ',I 
ITMzI 
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DO 222 KN-1,NMI 
CA-KN 
CB-CA«NMI 
READ(14,555) COEF(CA,ITM),COEF(CB,ITM),DUM 
WRITE(15,*) KN,COEF(CA,ITM),COEF(CB,ITM),DUM 


222 CONTINUE 


ITEZI-NA 
DO 333 KNz1,NMI 
CA=KN 
CB=CA+NMI1 
READ(14,555) COEF(CA,ITE),COEF(CB,ITE),DUM 
WRITE(15,*) KN,COEF(CA,ITE),COEF(CB,ITE), DUM 


333 CONTINUE 
665 CONTINUE 


€ 


Oi OO 


e 
C 


Skipping Field Generation for m=0 if +/- Z-Axis Incidence 
IF((M.EQ.0).AND.(IM1.EQ.1)) GO TO 777 


WRITE(*,*) 'Generate Fields Across Array' 


DX=RCORE/RPTS 
DY=DX 


DO 55 J=1,NPTS 
DO 33 L=1,NPTS 


Determine local planar array coordinates 


JI=FLOAT(J-1)-RPTS 
L1=FLOAT(L-1)-RPTS 
XP=J1*DX 
YP=L1*DY 

ZP=0.0 


If Array Tilted Then Transform Coordinates 
From local Planar array coordinates 
to global Core coordinates.... 


IF(IFLG.EQ.1) THEN 
XC=A(1,1)*XP+A(1,2)*YP+A(1,3)*ZP 
YC=A(2,1)*XP+A(2,2)*YP+A(2,3)*ZP 
ZC=A(3,1)*XP+A(3,2)*YP+A(3,3)*ZP 

ELSE 
XC-XP 
YC-YP 
ZCRZP 

ENDIF 


Determine if R <= RCORE, Else Skip Field Computation 


RSQ=XC*XC+YC*YC+ZC*ZC 
R=SQRT(RSQ) 
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IF (R.GT.RCORE) GO TO 32 

C Determine spherical Bessel function scaled radius argument 
RK=KR*R 

e Determine Complex Spherical Bessel function and Derivative 


NJN=M+NM1+MF 
CALL CSBSL(RK,NJN,CJ,DCJ) 


C Compute Core Coordinates and Functions 


THETA-ACOS(ZC/R) 
PHI=ARCTAN2(YC,XC) 
C=COS(THETA) 
S=SIN(THETA) 
CMP=EM*COS(M*PHD 
SMP-EM*SIN(M*PHI) 


Determine Legendre Polynomial and Derivative 


NLP=NM1+MF 
CALL LPAD(M,NLP,C,P,DP) 


Consider only the incident angle of interest 


I=NAI 
ITM=I 
ITE=ITM+NA 


Initialize additive KN-sum fields to zero 


EP1=(0.,0.) 
EP2=(0.,0.) 
ET 1=(0.,0.) 
ET2=(0.,0.) 
ER1=(0.,0.) 
ER2=(0.,0.) 


. . 6 ee ò a. òo 6 . ò% @ v 


DO 22 KN=1,NM1 
N=KN+M+MF-1 
NJ=N+1 
NL=KN+MF 
CA=KN 
CB=CA+NM1 


ET 1=ET1+COEF(CA,ITM)*M*CI(NJ)*P(NL) 
+COEF(CB,ITM)*DCIJ(NJ)*DP(NL) 
EP1=EP1+COEF(CA,ITM)*CJ(NJ)*DP(NL) 
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1 +COEF(CB ITM)*M*DCI(NJ)*P(NL) 
ER1=ER1+COEF(CB,ITM)*N*(N+1)*CJ(NJ)*P(NL) 


ET2-ET24COEF(CA,ITE)*M*CJ(NJ)*P(NL) 


1 +COEF(CB,ITE)*DCJ(NJ)*DP(NL) 
EP2=EP2+COEF(CA,ITE)*CJ(NJ)*DP(NL) 
1 +COEF(CB ITE)*M*DCI(ND*P(NL) 


ER2=ER2+COEF(CB,ITE)*N*(N+1)*CJ(NJ)*P(NL) 


. . òo ò òo @ @ @ ọọ òo ò ò ò% ò . 0 e 


ETH(J,L,1)=ETH(J,L,1)-JAY*SRU*CMP*ET1/R 
EPHI(J,L,1)=EPHI( L,1)+JAY*SRU*SMP*EP1/R 
ERAD(,L,1)=ERAD(J,L,1)-JAY*S*CMP*ER1/(SRE*RSQ) 


ETH(,L,2)=ETH(,L,2)}+SRU*SMP*ET2/R 
EPHI(J,L,2)=EPHI(J,L,2)+SRU*CMP*EP2/R 
ERAD(J,L,2)=ERAD(J,L,2)+S*SMP*ER2/(SRE*RSQ) 


32 CONTINUE 
33 CONTINUE 


C= End of loops across antenna plane 


777 CONTINUE 
€ 
C** End of modal loop for field analysis 


ee SRR KK HK le al ol 


C Determine magnitude and phase of Computed field in the direction 
C ofthe theoretical field vector, and the magnitude of the error. 
C 


SINA=SIN(ALPHA(NAI)) 
COSA-COS(ALPHA(NAD) 


C Loop through array to determine field components 


DO 634 J=1,NPTS 
DO 633 L=1,NPTS 


E Determine local planar array coordinates 
C 
J1=FLOAT(J-1)-RPTS 
L1=FLOAT(L-1)-RPTS 
=J1*DX 
YP=L1*DY 
ZP=0.0 
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If Array Tilted Then Transform Coordinates 
From Array to Core .... 


IF(IFLG.EQ.1) THEN 
XC=A(1,1)*XP+A(1,2)* YP+A(1,3)*ZP 
YC=A(2,1)*XP+A(2,2)* YP+A(2,3)*ZP 
ZC=A(3,1)*XP+A(3,2)*YP+A(3,3)*ZP 

ELSE 
XC-XP 
YC-YP 
ZC-ZP 

ENDIF 
RSQ-XC*XC«YC*YC4ZC*ZC 
R=SQRT(RSQ) 


Determine if R <= RCORE, Else Zero Fields 
IF (R.LE.RCORE) THEN 
Compute Core Coordinates and Functions 


THETA=ACOS(ZC/R) 
PHI=ARCTAN2(YC,XC) 
CT=COS(THETA) 
ST=SIN(THETA) 
CP=COS(PHI) 
SP-SIN(PHI) 


NOTE: R IS GIVEN IN TERMS OF KO*R 
== KO IS INCLUDED IN R, DEFINE KO=1.0 


KO=1.0 
TM and TE Plane waves are of unit magnitude 


EOTE=1.0 
EOTM=1.0 


Complex incident field is exponential(-j*ko*p.r) 


use field direction for incident angle 180 to 90 degrees 
EXP1= CEXP( -JAY*KO*( XC*SINA+ZC*COSA ) ) 


Compute plane wave field unit vectors and components 


E TM = cos(alpha)Ux - sin(alpha)Uz 


E TMT 
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EI(1) = COSA*CT*CP - SINA*(-ST) 
E(1) = EOTM*EXP1*EI(1) 

ETMP 
El(2) = COSA*(-SP) 
E(2) = EOTM*EXP1*EI(2) 

ETMR 
EI(3) = COSA*ST*CP - SINA*(CT) 
E(3) = EOTM*EXP1*EI(3) 

E TM TOT 

E(4)=SQRT(CABS(E(1))**2.0+CABS(E(2))**2.0+CABS(E(3))**2.0) 


E TE = Uy 


E TET 
EI(5) = CT*SP 

E(5) = EOTE*EXPI*EI(5) 

ETEP 

El(6) = CP 

E(6) - EOTE*EXPI*EI(6) 

E TER 

El(7) = ST*SP 

E(7) - EOTE*EXPI*EI(7) 
E TE TOT 

E(8)-SQRT(CABS(E(5))**2.0- CABS(E(6))**2.04-CABS(E(7))**2.0) 


Ak ol ok ok ak ok ak ak ak 


(6) 6) 


(sc AO Cs). 


Find magnitude and phase of Computed field in the direction 
of the theoretical field vector 
E . ei /Ei = (ethrETH+ephi*EPHI+erad*ERAD )/EO*EXPI 


EDOTI=EI(1)*ETH(J,L,1)+EI(2)*EPHI(,L,1)+EI(3)*ERAD(J,L,1) 
EDOT( L,1)=EDOT1/(EOTM*EXP1) 


EDOT2=EI(5)*ETH(J,L,2)+EI(6)*EPHI(J,L,2)+EI(7)*ERAD(JL,2) 
EDOT(J,L,2)=EDOT2/(EOTE*EXP1) 


Find rel. magnitude of Computed field perpendicular to the direction 
of the theoretical field vector, the error component in the 
Computed field. 

Rel error 2| E - (E.ei) eil / EO 


= sqrt{ [Eth - (E.ei)*ethl**2 4 [Erad - (E.ei)*eradl**2 
+ IEphi - (E.ei)*ephil**2 )/ EO 


Field error magnitude is computed 

EDTER 1= CABS(ETH(IL,1) -EDOTI*EI(1))* *2.0 
1 +CABS(EPHIG,L,1)-EDOT1*EI(2))**2.0 

1 +CABS(ERADG,L,1)-EDOT1*EI(3))**2.0 
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C 
C 


EDTERR(J,L,1)ÆSQRT(EDTER1)(ÆOTM) 


EDTER2= CABS(ETH(,L,2) -EDOT2*EI(5))**2.0 


+ CABS(EPHI(J,L,2)-EDOT2*EI(6))**2.0 
+ CABS(ERAD(J,L,2)-EDOT2*EI(7))**2.0 


EDTERR(,L,2)=SQRT(EDTER2)/(EOTE) 


ELSE 


R > RCORE , Set field components to zero 


EDOT(J,L, 1)=0.0 
EDOT(J,L,2)=0.0 
EDTERR(J,L,1)=0.0 
EDTERR(,L,2)-0.0 


ENDIF 


633 CONTINUE 
634 CONTINUE 


CEAESERE SS * 


WRITE(*,*) ‘SPHERICAL FIELDS COMPUTED’ 
WRITE(*,*) 


C Optional Listing of Coordinates and E-Fields 


C 


WRITE(*,*) 'Check Values of Coordinates and E-Field ? (Y/N):' 
READ(*,101) YN 


IF((Y N.EQ.'Y.OR.(YN.EQ.y)) THEN 


757 CONTINUE 


WRITE(*,*) ‘Enter Array Triple Index to Check (J,L,D' 
WRITE(*,*) ' J= x-point 1,2,..,,,NPTS 

WRITE(*,*) ' L= y-point 1,2,..,,,NPTS 
WRITE(*,*) ' I= 1 for TM, 2 for TE 

WRITE(*,*) 

READ(*,*) J,L,I 

EPM=CABS(EPHI(J,L,D) 

ETM=CABS(ETH(,L,D) 

ERM=CABS(ERAD(,L,D) 
ETOT=SQRT(EPM*EPM+ETM*ETM+ERM*ERM) 


Determine local planar array coordinates 


J1=FLOAT(-1)-RPTS 
L1=FLOAT(L-1)-RPTS 
XP=J1*DX 
YP=L1*DY 

ZP=0.0 


If Array Tilted Then Transform Coordinates 
From Array to Core .... 
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IF(IFLG.EQ.1) THEN 
XC=A(1,1)*XP+A(1,2)* YP+A(1,3)*ZP 
YC=A(2,1)*XP+A(2,2)* YP+A(2,3)*ZP 
ZC=A(3,1)*XP+A(3,2)* YP+A(3,3)*ZP 

ELSE 
XC=XP 
YC=YP 
ZC=ZP 

ENDIF 

WRITE(*,*) "XP, YP,ZP:' XP, YP,ZP 

WRITE(*,*) XC,YC,ZC:' XC,YC,ZC 

WRITE(*,*) 

WRITE(*,*) —E-RADI =' ERM 

WRITE(*,*) 'IE-THI ='ETM 

WRITE(*,*) ‘TE-PHI! =',EPM 

WRITE(*,*) "E-TOT! ='ETOT 

WRITE(*,*) 

WRITE(*,*) Look at Another Point ? (Y/N):' 

READ(*,101) YN 

IF((YN.EQ.'Y').OR.(YN.EQ.'y’)) go to 757 


ENDIF 


INEO oH Oo o goce ecce sek eoe sje eee eoe joe toe oe oec oc oe oe see see esos je e esee esce e se eee eie ie 


C Field representation selection menu 
ad ak ak ak ak ak ak ak ak ak ak ak ok ofc oc aic ac akc aic ac oic kc oc ac ok aic ac alc ac aic ac aic aic ajo 


C 
C 


700 WRITE(*,701) 


E 


* 


ak ak ak ak ak ak ak ak ak ak ak ak ak ak ak ak ak ak ak ak ak ak ak ak ak afe ake ak ak afc ak ad a ae a a ak ad a ak ak ak ak ake ak ae a ad ad ad al ade a ad ad ad ak a ad ade ale ale a a 


701 format(////////////5X  CORFLD OUTPUT GRAPHICS PRESENTATION'/8X, 
I'CORE FIELD REPRESENTATION SELECTION MENU //5X, 

2 Please select the type of field representation of interest 

BOX. 

4'CORE will display the following field representations: //8X, 


C 


51. 
6'2. 
po. 
844. 
95. 
1'6. 
Du 
3'8. 


4'9. 
51 


TM INCIDENT FIELD, E-THETA'/8X, 

TM INCIDENT FIELD, E-PHI '/8X, 

TM INCIDENT FIELD, E-RADIAL '/8X, 

TM INCIDENT FIELD, E-TOTAL FIELD MAGNITUDE '//8X, 
TE INCIDENT FIELD, E-THETA'/8X, 

TE INCIDENT FIELD, E-PHI '/8X, 

TE INCIDENT FIELD, E-RADIAL '/8X, 

TE INCIDENT FIELD, E-TOTAL FIELD MAGNITUDE '//8X, 
COMPARE COMPUTED FIELD TO INCIDENT FIELD'//8X, 

0. CHANGE ASPECT RATIO //8X, 


60. FINISHED WITH THIS ANGLE OF INCIDENCE //5X, 
6'Indicate your selection entering "1","2”,...,"10",or "0" 
T,5X) 


* 


ak ak ak ak ak ak ak ak ak ak ak ak ak ak ak ak ak ak ak ak e a ak 3k ai al ak al ak ad ai al al 


READ(* *) SELECT 
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Ce O AOS 


Check response to menu 


IF (SELECT.EQ.1) THEN 


Ae ale ade ole ale ai al ol ade a al al al ale ad ai al a a o a a a a a ae a al a al a a ak al al al a al ali a al a al ei ali il al oi oc oic ofc a a a al al ok ok 


DISPLAY TM INCIDENT FIELD, E-THETA 


CHECK FOR DESIRED REPRESENTATION OF FIELD, MAGNITUDE/PHASE 
WRITE(*,715) 


Me oe al ade al c oic oc al ade ale al ale al al ad a ade ade al a ade ai ae a al afe A ai ae ai afe ai al ai al ade a a afe ae al ae a al al ae al al a al a a a al a ai a a Ik a Ik Ik k k 


715 format///////1///5X,CORE FIELD PRESENTATIONS'/5X, 


E 


C 


1'You have selected the following field element //5X, 

21. TMINCIDENT FIELD, E-THETA'//8X, 

3'This field can be plotted in PHASE or MAGNITUDE, //5X, 
4'Please select the type of field representation of interest 
5'/8X, 

61. PHASE'/8X, 

72. MAGNITUDE /8X, 

8'Indicate your selection entering "1" or "2" 


9'J/5X) 


Me ade al al oe ole ade ale ae ade kc ok ale kc je oc oc je ale a ale a al oc oco ale ale al ae a al 3k al 


READ(*,*) SELCTI 
ASSIGN AS INDICATED BY RESPONSE TO MENU 


DO 712 J=1,NPTS 
DO 711 L=1,NPTS 
CPLOT=ETH(_L,1) 
IF (SELCT1.EQ.1) THEN 
PLOT(J,L)=(ARCTAN2(IMAG(CPLOT),REAL(CPLOT)))/DTR 
ELSE 
PLOT(J,L)=CABS(CPLOT) 
ENDIF 


711 CONTINUE 
712 CONTINUE 


C 


O 


LO CO 


IF (SELCT1.EQ.1) THEN 
TITLE2 2' TM INCIDENT FIELD, E-THETA PHASE' 
ELSE 


TITLE2 =' TM INCIDENT FIELD, E-THETA MAGNITUDE' 
ENDIF 


ak ak ak k ale ade ak ok ok ale al al ake ale ad ad afe ad afe ad afe ade a ok ae al ae af ok ie af ae al ae a al afe ahe afe afe afe afe afe afe afe afe afe afe Afe afe afe ae afe afe ai a ad ai ai a ai a k a ae 2c 


ELSEIF (SELECT.EQ.2) THEN 


ak ai oj aic oc oc oc oc oj iO iO RIOR iO IORI afe afe IR QR I RR IR I I ii II ROK ii ke ke >k 


DISPLAY TM INCIDENT FIELD, E-PHI 


CHECK FOR DESIRED REPRESENTATION OF FIELD, MAGNITUDE/PHASE 
WRITE(*,725) 


ak ak ac sk al ade ade ade ae ale Ak al ae ae ale ae ale aie a ae aie fe fe afe afe al ai a ade ae ai a ai a a al ad afe af al al al ok al a ai al a Ie Ie Ie ae a a Ie ak Ik ak ok ok ak a a 2k 
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725 format(////////////SX,;CORE FIELD PRESENTATIONS '/5X, 
1"You have selected the following field element' J/5X, 
22. TM INCIDENT FIELD, E-PHIT'//8X, 
3'This field can be plotted in PHASE or MAGNITUDE,'//5X, 
4'Please select the type of field representation of interest 
5'/8X, 
61. PHASE'/8X, 
72. MAGNITUDE /8X, 
8 Indicate your selection entering "1" or "2" 
9° //SX) 


C de ale ade ade ale ade ade ade ak ak ak ak ak ak ak ak ak ak ak ak ak ak ak a ae ae e ae a ae ai e a a 
READ(* *) SELCTI 
C ASSIGN AS INDICATED BY RESPONSE TO MENU 


DO 722 J=1,NPTS 
DO 721 L=1,NPTS 
CPLOT=EPHI( L,1) 
IF (SELCT1.EQ.1) THEN 
PLOT(J,L)=(ARCTAN2(IMAG(CPLOT),REAL(CPLOT)))/DTR 
ELSE 
PLOT(J,L)=CABS(CPLOT) 
ENDIF 
721 CONTINUE 
722 CONTINUE 


C 
IF (SELCT1.EQ.1) THEN 
TITLE2 ='TM INCIDENT FIELD, E-PHI PHASE ' 
ELSE 
TITLE2 = 'TM INCIDENT FIELD, E-PHI MAGNITUDE’ 
ENDIF 
C 
C e ae ae ad e c oe c oe oe oe oe oe ke ke oe ode aic ok aic ok ok ok ok oke ofc ok ok oke ac oke kc oc ok ok ode fe fe ae fe k ak fe okc ok ok ak okc ok ok ke oc oc aic aic oc ac oc oc ok oc aic oe oo 
ELSEIF (SELECT.EQ.3) THEN 


ak kc ke oc kc al oe ok occ ac kc oc oc ok ok ok oe oc oc oc ok oe ok oe oe oe ok a a a e e ae a a KK KK KKK KK KKK KKK KKK KK KK KKK 


C 
C DISPLAY TM INCIDENT FIELD, E-RADIAL 
€ 
C 


CHECK FOR DESIRED REPRESENTATION OF FIELD, MAGNITUDE/PHASE 
WRITE(*,735) 
ak ke ke c ke ke e oe o he Ak ke ad ae ae ae ad a a ae a ae ale a a a ae a ae ad a a ai a ae a 3k a e al a ae ae a a a a e e e e ae a Ak e e e e a a a 3k e a k k 
735 format(////////////SX,,;CORE FIELD PRESENTATIONS ',/5SX, 
1'You have selected the following field element' J/5X, 
23. TM INCIDENT FIELD, E-RADIAL'//8X, 
3'This field can be plotted in PHASE or MAGNITUDE; '//5X, 
4'Please select the type of field representation of interest 
SAX, 
61. PHASE'/8X, 
72. MAGNITUDE /8X, 
8'Indicate your selection entering "1" or "2" 
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9° ISX) 


Qo aoaaa akak aaka a ak ak aka ak ae a ae aa loko ook Xo ook 
READ(*,*) SELCTI 
C ASSIGN AS INDICATED BY RESPONSE TO MENU 


DO 732 J=1,NPTS 
DO 731 L=1,NPTS 
CPLOT=ERAD(.L,1) 
IF (SELCT1.EQ.1) THEN 
PLOT(J,L)=(ARCTAN2(IMAG(CPLOT),REAL(CPLOT)))/DTR 
ELSE 
PLOT(J,L)=CABS(CPLOT) 
ENDIF 
731 CONTINUE 
732 CONTINUE 


C 
IF (SELCT1.EQ.1) THEN 
TITLE2 ='TM INCIDENT FIELD, E-RADIAL PHASE’ 
ELSE 
TITLE2 ='TM INCIDENT FIELD, E-RADIAL MAGNITUDE’ 
ENDIF 
C 
C ak ak ade ade ae de a al ad al ad ak ak ak a ai a ake ak ak a a a a a ae al a a a a a al ak ak ak ake ak ak ak ake ake ak ak ak ak ak ak ak ak ak ak ak ak ak a a a a a ad a ad A a a 
ELSEIF (SELECT.EQ.4) THEN 
e ak ak ak ak ak ak ak ak ak ak ak ak ak 2k ak a oc kc ok ok ok ok aic ae ok ok abe ak ak ak ak ak ak ak ak ak ac ok ok aic okc aic aic ac ok ok aic aic ok aic oc aic ac ac aic aic aic ac abc aic ac adc oc xe eo 


C DISPLAY TM INCIDENT FIELD, E-TOTAL FIELD MAGNITUDE 
C ASSIGN AS INDICATED BY RESPONSE TO MENU 


DO 742 J=1,NPTS 

DO 741 L=1,NPTS 

EPM-CABS(EPHI(,L,1)) 

ETM-CABS(ETH(,L,1)) 

ERM-CABS(ERAD(,L,1)) 

PLOT(J,L)=SQRT(EPM*EPM+ETM*ETM+ERM*ERM) 
741 CONTINUE 
742 CONTINUE 


TITLE2 = "TM INCIDENT FIELD, E-TOTAL FIELD MAGNITUDE ' 


C ak ak akk ak ak ak ak ak akk ak ak a ad ai o akk ode ode ae ade al a a e a ak ake ake ak ak ak ake ak dd ai o a ai a a a a 2k a ake ake afe ake ak ak ak a ak 2k ak a a k ak k a a a a 


ELSEIF (SELECT.EQ.5) THEN 


ade ade ade ode ade ad ak ak ak ade ae ak akk akk akk ak akk ak ake ak ak afe ak ak afe ake e a a ake afe ake afe afe ok abc ac ae ac ac aic ac aic ac ac ak ac ak e a ok ok ok a ok ak ok a ad a ad Ak k k k 


DISPLAY TE INCIDENT FIELD, E-THETA 


CHECK FOR DESIRED REPRESENTATION OF FIELD, MAGNITUDE/PHASE 
WRITE(*,755) 


ak ak ak ak ak ak ak ak ak ak o oe a oe de oe a oe e ad al ak ak ad okc al ode ad al a a a al 2k okc oc ac ok 2k a a a oc oc ok ak ak a 2k ak ak ak al ak al ok oc ac oc oco oe k k 


uU C) uu 
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755 format(J///////////SX, CORE FIELD PRESENTATIONS'/5X, 
l'You have selected the following field element’ //5X, 
2'5. TE INCIDENT FIELD, E-THETA'//8X, 
3'This field can be plotted in PHASE or MAGNITUDE; //5X, 
4'Please select the type of field representation of interest 
5'/8X, 
61. PHASE'/8X, 
72. MAGNITUDE /8X, 
8 Indicate your selection entering "1" or "2" 


9'//5X) 


ak ak akk ak ade a a a a a a a a a a c c aic ode ac je ac aic ac aic ac aic alic ak aic ak ak ak ak 2k 


READ(*,*) SELCTI 


C ASSIGN AS INDICATED BY RESPONSE TO MENU 


751 


DO 752 J=1,NPTS 
DO 751 L=1,NPTS 
CPLOT=ETH( L.2) 
IF (SELCT1.EQ.1) THEN 
PLOT(,L)=(ARCTAN2(IMAG(CPLOT),REAL(CPLOT))/DTR 
ELSE 
PLOT(J,L)=CABS(CPLOT) 
ENDIF 
CONTINUE 


752 CONTINUE 


C 


C 


IF (SELCT1.EQ.1) THEN 


TITLE2 = ' TE INCIDENT FIELD, E-THETA PHASE' 
ELSE 


TITLE2 =' TE INCIDENT FIELD, E-THETA MAGNITUDE' 
ENDIF 


ak ak ak ak ak ak ak ak ak ak ak ak ak ak ak oic ac ak ake ak ake ak ake ake ak ake ak ak ak ake akc ak ak ak ak akc ak ak ak ak ak ak ake ak 3k ak ak ak ak ak ak 3K ak ak ak 3K ak ak ak a AR AR a a 


ELSEIF (SELECT.EQ.6) THEN 


ACA C 


ak ak ak ak ade a a akc ak ak ak a a ak ak a a c ak ake ake ak aic aic aic aic ac kc aic aic aic aic ac oc aic aic ae aic ac aic alc ac ac ac ac aic alc ac ac ac afe ac alc ac ac alc ac ac ac alc ac aic oc oo o9 


DISPLAY TE INCIDENT FIELD, E-PHI 


CHECK FOR DESIRED REPRESENTATION OF FIELD, MAGNITUDE/PHASE 


WRITE(*,765) 


ak ak ak ak ak ak ak ak ak ak ak ak ak ak ak aic ak ak ak ak aic aic aic alc c aic aic ac alc oic aic ac ac ok ak ak ak ak ak ak ak ak ak ak ak ak ak akk ak ak ak ak ak ak ak ak ak ak ak ak ak ak ak 3k 3k 2k 


765 format(///////[////5X,'| CORE FIELD PRESENTATIONS'/5X, 
l'You have selected the following field element //5X, 
2'6. TE INCIDENT FIELD, E-PHIT'//8X, 
3'This field can be plotted in PHASE or MAGNITUDE, //5X, 
4'Please select the type of field representation of interest 
5'/8X, 
61. PHASE'/8X, 
72. MAGNITUDE /8X, 
8'Indicate your selection entering "1" or "2" 
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9'//SX) 


C ak ak kc ak ok ok ok ok ok ok ok ok ok ae ae ae a ok ok ok ae a ad ad a ok ok a a a ae al ae ad 
READ(*,*) SELCTI 
C ASSIGN AS INDICATED BY RESPONSE TO MENU 


DO 762 J=1,NPTS 
DO 761 L=1,NPTS 
CPLOT=EPHI(J,L,2) 
IF (SELCT1.EQ.1) THEN 
PLOT(J,L)=(ARCTAN2(IMAG(CPLOT),REAL(CPLOT)))/DTR 
ELSE 
PLOT(J,L)=CABS(CPLOT) 
ENDIF 
761 CONTINUE 
762 CONTINUE 


E 
IF (SELCT1.EQ.1) THEN 
TITLE2 ='TE INCIDENT FIELD, E-PHI PHASE’ 
ELSE 
TITLE2 ='TE INCIDENT FIELD, E-PHI MAGNITUDE ' 
ENDIF 
€ 
C ak ak ak ok ok ak ok ak ok ok ok ok ok ok ok ok ok ok ok ok ok ok ok ok ok ok ok ok ok ok ok okc ok ok ok ok ok ok ok ok ok ok akc ok ok ok ok ok ok ok ok ok ok ok ok ac ok ok ok ok ale ae a oc oco 
ELSEIF (SELECT.EQ.7) THEN 
C sk ak ak ok kc ok ok ok ok ok okc ok ok ok ok ok ok okc ok ok ok ok ok ok ok okc ok ok ok ok ok ok ok ok ok ok ok ok ok ok ok ok ok ok okc ok ok ok ok ok ok ok ok ok ok okc ok ok ok ok ok ok ok ok ad al 
C DISPLAY TE INCIDENT FIELD, E-RADIAL 
C | CHECK FOR DESIRED REPRESENTATION OF FIELD, MAGNITUDE/PHASE 
WRITE(*,775) 
e ak ak ok ok ok ok ok ok ok ok ok ok ok ok ok ok ok ok ok ok ok ok ok ok ok ok ok ok ok ok ok ok ok ok ok ok ok ok ok ok ok ok ok ok ok ok ok ok ok ok ok ok ok ok ok ok okcok ok ok ok ok ok ok oc ok 


TTS format(////////////5X ,' CORE FIELD PRESENTATIONS'/5X, 
l'You have selected the following field element //5X, 
27. TEINCIDENT FIELD, E-RADIAL'//8X, 
3'This field can be plotted in PHASE or MAGNITUDE, //5X, 
4'Please select the type of field representation of interest 
5'/8X, 
61. PHASE'/8X, 
72. MAGNITUDE /8X, 
8'Indicate your selection entering "1" or "2" 
9'//5 X) 


C ak ak kc ok ok ok ok ok ok ok ok ok ok ok ok ok ok ok ok ok ok ok ok ok ok ok okc ok ok ok ok ade ad ale 
READ(*,*) SELCTI 


C  ASSIGN AS INDICATED BY RESPONSE TO MENU 


DO 772 J=1,NPTS 
DO 771 L=1,NPTS 
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CPLOT=ERAD(J,L,2) 
IF (SELCT1.EQ.1) THEN 
PLOT(J,L)=(ARCTAN2(IMAG(CPLOT),REAL(CPLOT)))/DTR 


ELSE 
PLOT(,L)=CABS(CPLOT) 
ENDIF 
7711 CONTINUE 
772 CONTINUE 
E 
IF (SELCT1.EQ.1) THEN 
TITLE2 = ' TE INCIDENT FIELD, E-RADIAL PHASE' 
ELSE 
TITLE2 =' TE INCIDENT FIELD, E-RADIAL MAGNITUDE' 
ENDIF 
C 
Ç aaka aao oaaao aaao aooaa kkk kkk kokk 
ELSEIF (SELECT.EQ.8) THEN 
Q^ okkeokekookokek okookokokookok aoaaa oaaao aaoo kk 


C DISPLAY TE INCIDENT FIELD, E-TOTAL FIELD MAGNITUDE 


DO 782 J=1,NPTS 
DO 781 L=1,NPTS 
EPM=CABS(EPHI(J,L,2)) 
ETM=CABSŒTH(J,L,2)) 
ERM=CABS(ŒRAD(J,L,2)) 
PLOT(J,L)=SQRT(EPM*EPM+ETM*ETM+ERM*ERM) 
781 CONTINUE 
782 CONTINUE 
TITLE2 ='TE INCIDENT FIELD, E-TOTAL FIELD MAGNITUDE ' 


C 3k ak a 3k ak ak ak kc ok ok kc oc ok ke ok ok ok ok ok ak ak ak ak ak ak ak ak ok kc oc oc ok ak ak ak ok ke ok ale ad a ae al ak ak ak oc ok ke ok ae ad a oe ke ke ok ok ak ad oe ok oe ke x 


ELSEIF (SELECT.EQ.9) THEN 
C ak ak ak ak ak ak ak ak ak ak ak oc ok oc ok ae a ak a oe oc oc oc ok ok oc ok ak ak oc ok oc ok ok oko ok ak ok ak ok oe ke ook ok ok a ok ok ak ak ok oc ok ak ak ak ak ak ak ak ak ak ak ak ak ak ak 3k a 
C 
C Incident Field comparison representation selection sub-menu 
C ak ak ak ak c e oe ak ak ak ak ak ak ak ak ak ak ak ak ak ak ak ak ak ak ak ak ak ak ak ak ak ak ak ak ak ak ak ak ak ak ak ak ak ak ak ak ak ok ok ook ok ok ok ok ok ok oe o 


© 


800 WRITE(*,801) 
C ak ke ak c oe e e ke o ke oe oe e oe eo oe oe oe oe oe ke ok ok oc ok oe ok oc oe oc oe ok ak ok ok ak ak ak oe oe o ak ak 3k ak ak ok o ak ok ok ak oc ok ak oc oc oc oe oe oe oe oe oe 

801 format(////////////SX ,' CORE COMPARISON TO IDEAL FIELD'//5X, 
1'CORE FIELD REPRESENTATION SELECTION MENU 2'//5X, 
2'Please select the type of field representation of interest 
SU, 
4'CORE will display the following field representauons:' //8X, 
51. TMINCIDENT FIELD, Computed Field dot Ideal Field'/8X, 
62. TM INCIDENT FIELD,  Scaled Error component '///8X, 
73. TE INCIDENT FIELD, Computed Field dot Ideal Field'/8X, 
84. TE INCIDENT FIELD, Scaled Error component '///8X, 
9'Indicate your selection entering "1","2","3", or "4" 
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1'//5X) 


C e ale al ale k ak ak ak ak ak ak ak ak ak ak ak Ak ak ak ak ak ak ak ak ak ak ak 3k k k k k k k 


READ(*,*) SELCT2 


C Check response to menu 


A ak ak a o e 


Quo 


IF (SELCT2.EQ.1) THEN 
ak ak ade ade ade ak ake ake 3k ak ak 3k Ak ak a Ak ak ak ae a e de ak e a ak ak ak ak aic ak ak ak ak ak ak ak ak ak ic c ic ic Ak ak ic i aic ic e a a a ak ak ak e e a a k k k k k 
DISPLAY TM INCIDENT FIELD, Computed Field dot Ideal Field 


CHECK FOR DESIRED REPRESENTATION OF FIELD, MAGNITUDE/PHASE 
WRITE(*,815) 


de ade ade a ade al a a a a a a a ad a a a a e a a e a a e a a e a a e a a a de a id e e a a a a a a e a a HK HK HK KK KK KK KK EK KH 


815 format(////////////5 X CORE FIELD PRESENTATIONS'/5X, 


E 


e 


1"You have selected the following field element //5X, 

2'1. TM INCIDENT FIELD, Computed Field dot Ideal Field’ //8X, 
3'This field can be plotted in PHASE or MAGNITUDE; '//5X, 
4'Please select the type of field representation of interest 

5'/8X, 

61. PHASE'/8X, 

72. MAGNITUDE /8X, 

8'Indicate your selection entering "1" or "2" 

9' //5X) 


ak jc seb dc je je se e o e oc bee c e a jc ade oc a c ok oe ak jc 2k ake ak 2k 
READ(*,*) SELCTI 


ASSIGN AS INDICATED BY RESPONSE TO MENU 


DO 812 J=1,NPTS 
DO 811 L=1,NPTS 
CPLOT=EDOT(.L,1) 
IF (SELCT1.EQ.1) THEN 
PLOT(,L)=(ARCTAN2(IMAG(CPLOT),REAL(CPLOT)))/DTR 
ELSE 
PLOT(,L)=CABS(CPLOT) 
ENDIF 


811 CONTINUE 
812 CONTINUE 


C 


C 


IF (SELCT1.EQ.1) THEN 


TITLE2='TM INCIDENT FIELD, Computed Field dot Ideal Field, PHASE' 
ELSE 


TITLE2="TM INCIDENT FIELD, Computed Field dot Ideal Field, MAGN.’ 
ENDIF 


Mee e ad e c oc ak ak a ak ak ak ad oc o al a a a a a a a a al a a a a a oc oc oc ac oc oc ac ac okc oc e e a a ak a ak a ak k ak a a ak OK KK KE k 


ELSEIF (SELCT2.EQ.2) THEN 


ale ade ak a k ak ak ak lalala a a a e de ad a a e fl al al de ad a e al al de dl al al e all AAA 


148 


C DISPLAY TM INCIDENT FIELD, Scaled Error component 


DO 822 J=1,NPTS 
DO 821 L=1,NPTS 
PLOT(J,L)=EDTERRUJ,L,1) 
821 CONTINUE 
822 CONTINUE 
TITLE2 = 'TM INCIDENT FIELD, Scaled Error component ' 


ak ak ak ake de a a a al al al a ad al a ae a ai ai a a a a a ae ai a a a al ak a$ akk ake ake akk afk a a a a al akk ake akk a al akk al ae akc ac ake ake a ai ake ak a a ad ad a ak ak ak 
ELSEIF (SELCT2.EQ.3) THEN 
ak ak ak P SESS SES ES STE SS ESE SE ESE SES TES STE ST STS ST SS TCS TC SST TT SS TT TS SS TE EF 
DISPLAY TE INCIDENT FIELD, Computed Field dot Ideal Field 
CHECK FOR DESIRED REPRESENTATION OF FIELD, MAGNITUDE/PHASE 
WRITE(*,815) 
ak ak ak ak ak ak ak ake ak a ak ake ak ak ok okc ok ok ok al ae ad ae a a a a al ad a a a ad a ai al a a a ad a a a a a ali a ae a ae a a a ae a ai a a a a a a a ae a a 
835 format(/////1/////5SXCORE FIELD PRESENTATIONS'/5X, 
1"You have selected the following field element’,//5X, 
2'3. TE INCIDENT FIELD, Computed Field dot Ideal Field' //8X, 
3'This field can be plotted in PHASE or MAGNITUDE, ',//5X, 
4'Please select the type of field representation of interest 
5'J8X, 
61. PHASE /8X, 
72. MAGNITUDE /8X, 
8'Indicate your selection entering "1" or "2" 
9'//5X) 


c ak ak ak ak ak ak ak ak ak ak ak al ad ake ak al ad ake ak a al ad ai a a oic oco a al al ak 3k 


O AO 


READ(* *) SELCTI 
C ASSIGN AS INDICATED BY RESPONSE TO MENU 


DO 832 J=1,NPTS 
DO 831 L=1,NPTS 
CPLOT=EDOT(.L 2) 
IF (SELCT1.EQ.1) THEN 
PLOT(J,L)=(ARCTAN2(IMAG(CPLOT),REAL(CPLOT)))/DTR 
ELSE 
PLOTQ,L)-CABS(CPLOT) 
ENDIF 
831 CONTINUE 
832 CONTINUE 
C 
IF (SELCT1.EQ.1) THEN 
TITLE2='TE INCIDENT FIELD, Computed Field dot Ideal Field, PHASE' 
ELSE 
TITLE2='TE INCIDENT FIELD, Computed Field dot Ideal Field, MAGN.’ 
ENDIF 
E 


C 3k ak ak ak ak ak ak ab ab ok ak ad ak akk akk ak al ade ade ae ai ade al akk ak ak akk ake akk akk ae ad a ad a ade ad al ad kc kc okc ok ak ake ak akk okc ok okcok ok afk ak ak akk F akk ak ad e a a ad OK ok 
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ELSEIF (SELCT2.EQ.4) THEN 


C ak ak ak akak ak ak ak ak Ak ak ak ak ak ak ak ak ak ak ak ak ak Ak oc aie oc oc ac ac oc aie oc oc oic oc oc oc oc oic oic ac oic alc ac oic oic oc ac oc ok oie oc oc o ak ak ak ak ak ak ak ak ak ak ak ok 


C DISPLAY TE INCIDENT FIELD, Scaled Error component 


DO 842 J=1,NPTS 
DO 841 L=1,NPTS 
PLOT(J,L)=EDTERR(,L,2) 
841 CONTINUE 
842 CONTINUE 
TITLE2 = 'TE INCIDENT FIELD, Scaled Error component ' 


C ak al ode ade ade ade de ad al ale ale ok ok al ad ak de ad ad ad ak al ak a al ok ak ok ak ak ok ak ak ak ak ak ak ak ll ad ak ok ok ak ak ok ok ok a ad ok ode oe oe i oe oe oe oe oe a ok ok ok 
ELSE 
e ak ak ak ak ak ak ak ak ak oc oe ok oc oc oc oc oc ok oe oco ok oc oc ok oc oc ok oko Ak Ak Ak ak ak ak ak ak ak ak oe ok ok ak ak ak ok ook ak ak ak ak ak ak ak ak ak ak ak oc oc oc ok k 
C  ((SELCT2LT.0).OR.(SELCT2.GT.4) ) 
GO TO 800 
ENDIF 
E 


C ak ak ak Ak ak ak al al ak ak a ad ak ak ak ak Ak al ak ak ak ak Ak a Ak ak ak ak oc ok ak ak ak ak ak ak ak ak ak ak oc ok ak oc ok ak ak Ak ak Ak ak ak ak ak ak ak ak Ak ak ak ak Ak oc oe oe oe oe oe 


ELSEIF (SELECT.EQ.10) THEN 
e ak ak ak ak ak ok o ak ak oc oc oe oc oc oc oc ok ak ak ak Ak Ak o oc o ak ak Ak ak ak ak ak ak ak oi ak ak ak ak ak ak ak ak ak ak ak oi oc oc oc oc oc ok oc ok ac oc oc oc oc oc ok o oc o ok 
C ASPECT RATIO 
791 write(*,792) 
792 FORMAT(/SX,"***###*##* ASPECT RATIO *# #48884 #! fJ x 
l'IN ORDER TO AVOID DISTORTION ON THE SCREEN, IT IS NECESSARY’ /5X, 
2'TO DEFINE THE NUMBER OF ROWS PER INCH AND THE NUMBER OF'/5X, 
3'COLUMNS PER INCH. THESE DEFINITIONS CHANGE AS THE'/5X, 
4'RESOLUTION OF THE SCREEN, AS WELL AS THE TYPE OF SCREEN'/5X, 
S'CHANGE. //5X, 


C "THIS PROGRAM WAS WRITTEN IN MODE 16 (EGA) USING'/5X, 
C 6AGB-1 VIDEO BOARD AND A NEC MULTISYNC MONITOR. THE'/5X, 
C 7'OPTIMUM ASPECT RATIO FOR THIS CONFIGURATION WAS .65'///5X, 
8'PLEASE INPUT ASPECT RATIO: //5X) 
read(*,*) aspect 
goto 700 
C ak ak ak ak ak ak ak al ak ak ak ak akak ak ak ak ak ak ak ak Ak a ak ak ak ak ak ak Ak ak ak ok ok ak ak ak ak ak ak ak ak ak ak ak ak ak e al ak Ak ak ak ak ak ak Ak ak ak ak ak ak ak ok 3 k 
ELSEIF (SELECT.EQ.0) THEN 
C ak ak ak ok ak ak ak ak al Ak ak ak a al oc ok oe ok ak ak ak ak ad ak a ak a ak ak k ak al a ak ak ak ak Ak a ak ak ak Ak ak Ak Ak Ak ak ok ak ak ak Ak ak oe oe oc oko oko k k k 
C FINISHED WITH THIS ANGLE OF INCIDENCE 
goto 999 
C ak ak ak oc oc ok ak ak ak ak ak ak ak ak ak ak ak ak ak al al ak ok Ak ak ak ak ak Fk ak ak ak ak ak ak ak ok akak ak ak ak ak Fk ak Ak ak ak Ak ai ak ak ak ak ak ak ak ak ak ak ok ak ak ake akc ok 
ELSE 
C ak ak ak ak ak ak ak ak ak ak ak ak ak ak ak ak ak Ak ak ak ak Ak ak ak ai ak ak ak ak ak ak ak ak oc okc ak ak ak ak ak Ak ak ak ak ak ak ak ak ak ak ak ak ak ak ak Ak ak ak ak ak ak ak F akc ok 


C  ((SELECT.LT.0).OR.(SELECT.GT.9) ) 
GO TO 700 
ENDIF 
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CC 2G GIG GIGI aaka aaa a akak ak ak ak ak akak ak ak ak ak akak ak ak ak ak ak ak ak ak ake ak ake ak ak ak ak ak a ake ake ak ake ak ak ak ak 2k ak 


C 


WRITE (*,*) ' CALL PLOTTING ROUTINE' 
CALL PLOT3D(PLOT,50,50,ASPECT,GRFLAB, 

+ TITLE 1,TITLE2,HDR1,HDR2,ANGINC(NAD,DTH) 
GO TO 700 


C End of graphic output routine 


[EN AER AN JN ic ic sje xci xicsjexje exc sexe jeokc sje le aleae ac cae eje cle ok e je ic o s ee ok ARA EAE AHH 


999 CONTINUE 
C 
C Program completed 
C 
C I/O FORMAT STATEMENTS 
C 


10 FORMATV(//IlIIIIII[I[[/T1X, * **** WELCOME TO CORFLD *****' //7X , 
1"THIS PROGRAM COMPUTES INTERIOR EM FIELD WITHIN PENETRABLE /7X, 
2 BODIES OF REVOLUTION USING THE FIELD COMPONENT INFORMATION /7X, 
3'OUTPUT BY THE PROGRAM EMCAD.'/7X, 
4'PLEASE PRESS ANY KEY TO CONTINUE.',/) 

14 FORMATUQ//IMIMIMIIIIIIT TX, 
I'THE OUTPUT DATA FILE IS THE FILE CONTAINING THE OUTPUT /7X, 
2RESULTS OF ALL NUMERICAL CALCULATIONS CONDUCTED BY EMCAD.'/7X, 
3'THE FORMAT FOR THIS INPUT IS FILENAME ONLY. NO EXTENSION IS'/7X, 
4'REQUIRED OR DESIRED. EMCAD AUTOMATICALLY APPENDS AN '/7X, 
S'EXTENSION OF .INT TO YOUR FILENAME. //7X, 
6 PLEASE ENTER THE FILENAME OF THE OUTPUT DATA FILE. //7X) 

23 FORMAT(/7X, ENTER INC FLD ANGLE (DEG) FOR # ',13,//7X) 

25 FORMAT(//////7X,, PLEASE PRESS ANY KEY TO CONTINUE. ////) 

100 FORMAT(/,"* 2% 888 #88 8 CORE OUTPUT DATA ***soockooerororor 
arkit) 

101 FORMAT(A) 

102 FORMAT(I5) 

103 FORMAT(4(E14.6)) 

104 FORMAT( *** PROGRAM ABORTED BECAUSE -------- TER) 

107 FORMAT(/7X,15,22X,(,1PE11.3,2X,1PE11.3,))) 

109 FORMAT(/7X, INCIDENT FIELD ANGLES'//9X,'N',.10X, ALPHA(N)) 

110 FORMAT(/X,I3,7X,F5.0, DEG) 

111 FORMAT(/7X, COMPLEX Er(n) AND Ur(n)) 

112 FORMAT(//7X, ENTER INTERIOR DATA FILE (D:FILENAME.EXTENSION): '’) 

115 FORMAT(//7X, ENTER CAPTION OR LABEL ) 

124 FORMAT(/7X, INDEXING PROGRAM THROUGH VALUES OF M?) 

125 FORMAT(7X,'M-LOOP .... M = 5) 

135 FORMAT(7X,EX M-LOOP) 

136 FORMAT(//7X,'********* CORE ANALYSIS COMPLETED **************') 

204 FORMAT(J/7X , *oeieooooookook M o 1 S, olookokolokiokelelorior n 

300 FORMAT(7X,15,5E12.3) 
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500 FORMAT(3X,15,A64) 

503 FORMAT(3X,E14.6,A64) 

505 FORMAT(3X,15,2X,15,A64) 

507 FORMAT(3X,1PE11.3,2X,1PE11.3,2X,A64) 
510 FORMAT(3X,F5.0,A64) 

511 FORMAT(3X,F5.0) 

555 FORMAT(3X,E12.6,3X,E12.6,4X,E12.6,3X,E12.6,3X,A64) 
566 FORMAT(3X,A10,E12.6,3X,E12.6) 

567 FORMAT(3X,A10,E12.6) 

576 FORMAT(3X,A5,15,A2,E12.6,3X,E12.6) 
577 FORMAT(3X,A5,E12.6) 

578 FORMAT(3X,E12.6) 


99 CONTINUE 
C 
C Program termination 
C 
WRITE(*,136) 


STOP 
END 


ak al ak ak ak al ae ade a ak oc i oc oe be ok ok ake ake oc oe oe oc oe oc oe oc oc oc oe oc ok oc oc oe oe oc ok oc oc oe oe oc oe oc oc oc ok oc oe oc oc oc oc oc oe oe oe oe oe eoe 


SUBROUTINE CSBSL (Z,N1,CJ,DCJ) 
COMPLEX Z,Z2,F1,F2,FSTO,ALPHA,CJ(*),DCJ(*) 


This routine computes Riccati spherical bessel functions and 
derivatives for complex Z using downward recurrence as in 
Miller's algorithm, but with highly accurate starting values 
using truncated Tatlor's series. 


Returns CJ(n)=J[n-1] (kr) for nz1...N1 
DCJ(n)-J [n-1] (kr) 


DEFINING CONSTANTS 


SA C) 


ZM2=(CABS(Z))**2 
NMX=N1+ZM2 
Z2=(Z*Z)/2.0 
D1=2.*NMX+3. 
D2=D1*(2.*NMX+5.) 
D3=D2*(2.*NMX+7.) 
D4=D3*(2.*NMX+9.) 


Using Truncated Taylor's series for relative upper starting 
values of J(NMX) and J(NMX+1) with NMX .GT. ABS(Z*Z). 


GOA 


F1=1.-Z2/D1+Z2*Z2/(2.*D2)-Z2*Z2*Z2/(6.*D3) 


t52 


F2=Z*(1./D1-Z2/D2+Z2*Z2/(2.*D3)-Z2*Z2*Z2/(6.*D4)) 


C Performing downward recurrence with highest two orders in CJ(N1). 
C 
M=NMX 
11 FSTO=F1 
F1=(2.*M+1.)*F1/Z-F2 
F2=FSTO 
IF (CABS(F1).LT.1.0E20) GO TO 1 
F1=F1*1.0E-20 
F2=F2*1.0E-20 
FSTO=FSTO*1.0E-20 
1 CONTINUE 
M=M-1 
IF (M+2.GT.N1) GO TO 11 
CJ(N1)=F2 
CJ(N1-1)=F1 
NO=N1-2 


C Continuing downward recurrence to CJ(1) 
C 


DO 22 K=1,N0 
J=N1-K-1 
22 CJ(J)=(2.*J+1.)*CI(I+1)/Z-CJ (J+2) 


C Normalizing entire array wrt actual JO(Z). 
€ 
ALPHA=CSIN(Z)/CJ(1) 
DO 33 K=1,N1 
33 CJ(K)=ALPHA*CJ(K) 


C Computing derivatives 
C 
DCJ(1)=CI(1)/Z-CI(2) 
DO 44 K=2,N1 
44 DCK(K)-CJ(K-1)-(K-1)*CY(K)/Z 
RETURN 
END 


Ad ak ak al ak ak ak ak ade ad ak ak ak ak ak al ad ak ak ak ad ak ak ak ak ak akk ak ak akc ak ak ak ak ak 


SUBROUTINE LPAD (M,N1,C,P,DP) 


This routine computes Normalized Legendre polynomials and 
derivatives. 


For m >= 1; P(n) = P(m,n+m-1)/sin(THETA) 

DP(n) = dP(m,n+m-1)/dTHETA for nz1,2,... 
For m = 0; P(n) = P(0,n)/sin( THETA) 

DP(n) 2 dP(0,n) / ATHETA for nz1,2,... 


OHO Oi ©) 


REAL P(*),DP(*) 
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KSP=N1-1 
Q=1.0-C*C 
S=SQRT(Q) 
IF(M.GT.0) GO TO 1 
P(1)=1.0 

P(2)=C 

DP(1)=0.0 

DP(2)=-S 

Q=S 

GOTO 3 

1 IF (M.EQ.1) GO TO 5 
P(1)=-((-S)**(M-1)) 
DP(1)=-M*C*(-S)**(M-1) 
GOTO7 

5 P(1)=-1.0 
DP(1)=-M*C 

7 DO 2 K=1,M 
P(1)=P(1)*(2*K-1) 

2 DP(1)=DP(1)*(2*K-1) 
P(2)=P(1)*C*(2*M+1) 
DP(2)=(C*DP(1)-Q*P(1))*(2*M+1) 

3 DO 4 K=2,KSP 
N=K+M-1 
AK=M-N-1 
P(K+1)=((M+N)*P(K-1)-(2*N+1)*C*P(K))/AK 

4 DP(K+1)=((M+N)*DP(K-1)+(2*N+1)*(Q*P(K)-C*DP(K)))/AK 
RETURN 
END 


FUNCTION ARCTAN2(Y,X) 


Routine returns ATAN(Y/X) for all cases 


(5696) 


REAL ARCTAN2,Y,X PI 
PI=3.1415927 
IF(Y.NE.0.0) GO TO 11 
IF(X.EQ.0.0) GO TO 33 
IF(X.GT.0.0) ARCTAN2-0.0 
IF(X.LT.0.0) ARCTAN2-PI 
RETURN 

11 IF(X.NE.0.0) GO TO 22 
IF(Y.GT.0.0) ARCTAN2=0.5*PI 
IF(Y.LT.0.0) ARCTAN2-1.5*PI 
RETURN 

22 ARCTAN2=ATAN2(Y,X) 
RETURN 

33 ARCTAN2=0.0 
RETURN 
END 


E ak ak ak ae ade ad ak ak ak a ak a a ake ake ake akk ak ke oc okc okc okc okc ok ad ad ae ae ae a a ac okc alc okc okc ok aj ade ade ade ad okc ok a ae ae ad ad a ae aj ae ae ad aj aj aj ale ale a OK a 


C PLOT3D Plotting subroutine 
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REVISED: 22 NOV 1989 
This subroutine automatically scales and plots data as 3-D surfaces 


C 
C 
C 
€ 
C All subroutines beginning with 'Q' 

(t are from the Grafmatic plotting library 

C 

SUBROUTINE PLOT3D(Z,M,N,ASPECT,GRFLAB, 

+ TITLE 1, TITLE2,HDR1,HDR2,ALPHA,GAMMA) 

REAL X(50,50),Y(50,50),2(50,50), X MIN,XMAX,Y MIN, YMAX,ZMAX,ZMIN 
REAL XST, XFIN, Y ST, YFIN,ZST,ZFIN,ZFLOOR,ZAVG 

REAL PHI,THETA,ALPHA,GAMMA,STHETA,I1 J1,M2,N2,R 

INTEGER M,N 

CHARACTER*1 YES, GRAFPRNT 

CHARACTER*5 ASTRING,GSTRING,TSTRING,PSTRING,MIN,MAX,AVG 
CHARACTER* 12 name,label 

CHARACTER*64 TITLE1, TITLE2, TITLE3, TITLE4, GRFLAB, HDR1, HDR2 
CHARACTER*64 HCOPY 

DIMENSION IWORK1(640),IWORK2(640) 


C SETTING CONSTANTS 
€ 

BELL=CHAR(7) 

FF = CHAR(12) 


C Set constant for degree to radian conversion 
(C 

PI 2 3.14159265359 

DTR= PI/180. 


C Calculate # of points on array 
C 
NARRAY=PI*( (M/2)**2) 
M2=FLOAT(M-1)/2.0 
N2=FLOAT(N-1)/2.0 


IWORKI and IWORK2 should have a dimension equal to the number of 
pixels (or points on a pen plotter) in the x direction. 


O 2, 


NPTS=640 


C NPTS is the size of the IWORK1 and IWORK2 arrays 
C 
MODEP=16 


C Plot will be in graphics mode 16 
C 


C XYSET is a routine called to set up the initial x,y values (see below) 
C 
CALL XYSET(X,Y,M,N) 
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XST = X(1,1) 
YST = Y(1,1) 
ZST = Z(1,1) 


XFIN = XST 
YFIN = YST 
ZFIN = ZST 


ZMIN = 1000 
ZMAX = -1000 
ZAVG = 0.0 


PHI = 45.0 
THETA = 75.0 


DO 101=1,M 
DO 15J=1N 
IF (XJ) LT. XST) XST = XJ) 
IF (X(LJ) .GT. XFIN) XFIN = X(1J) 


IF (Y(1J) .LT. YST) YST - Y (1J) 
IF (Y(1J) .GT. YFIN) YFIN = Y(I,J) 


IF (Z(J) .LT. ZST) ZST - Z(1J) 
IF (Z(1J) .GT. ZFIN) ZFIN - Z(1JJ) 


C Check if point is on the antenna array 


11=FLOAT()-1.0 
J1=FLOAT(J)-1.0 
R= SQRT( ((11-M2)/M2)**2.0 + ((J1-N2)/N2)**2.0 ) 
IF (R.LT.1.0) THEN 
IF (Z(1J) .LT. ZMIN) ZMIN = Z(1,J) 
IF (Z(1J) .GT. ZMAX) ZMAX - Z(1J) 
ZAVG = ZAVG + Z(1,J)/(FLOAT(NARRAY)) 


ENDIF 
15 CONTINUE 
10 CONTINUE 


XFIN = XFIN + .2 
YFIN = YFIN + .2 


C The following parameters set the axes values, etc. (see documentation 
C for Q3DXAX, Q3DYAX, Q3DZAX) 
C 

XBEG=XST 

YBEG=YST 

XEND=XFIN 

YEND=YFIN 


ASTEP= 7 
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YSTEP= 2 
IF ((XFIN-XST).GE. 0.9) XSTEP =.5 
IF ((YFIN-YST).GE. 0.9) YSTEP = .5 


ZDIF=ZFIN-ZST 


IF (ZDIF.LE. 0.3) THEN 
ZSTEP = 0.05 
ZFIN = FLOAT(INT(10*ZFIN))/10+2*ZS TEP 
ZST = FLOAT(INT(10*ZST))/10-ZSTEP 
ZMINOR = 1 

ELSEIF (ZDIF.LE. 1.0) THEN 
ZSTEP = 0.25 
ZFIN = FLOAT(INT(ZFIN+0.5))+2*ZSTEP 
ZST = FLOAT(INT(ZST-0.5))-ZSTEP 
ZMINOR = 1 

ELSEIF (ZDIF.LE. 2.0) THEN 
ZSTEP = 0.5 
ZFIN = FLOAT(INT (ZFIN+0.5))+2*ZSTEP 
ZST = FLOAT(INT(ZST-0.5))-ZSTEP 
ZMINOR = 1 

ELSEIF (ZDIF.LE. 5.0) THEN 
ZSTEP = 1.0 
ZFIN = FLOAT(INT(ZFIN+0.5))+2*ZSTEP 
ZST = FLOAT(INT(ZST-0.5))-ZSTEP 
ZMINOR = 1 

ELSEIF (ZDIF.LE. 10.0) THEN 
STEP = 2.0 
ZFIN = FLOAT(INT(ZFIN+0.5))+2*ZSTEP 
ZST = FLOAT(INT(ZST-0.5))-ZSTEP 


ZMINOR = 1 
ELSE 
ZSTEP = 30.0 


ZFIN = INT(ZFIN+ZSTEP) 
ZST = INT(ZST-15.0) 
ZMINOR = 1 

ENDIF 

ZFLOOR=0.0 


C Clear the screen and enter input data 
E 
CALL QSMODE(2) 


C 3D-stick is a “wire-frame” figure with hidden lines removed. 

C 3D-fill is similar, except that the interior regions of each 

C "panel" making up the surface if filled in with color "klrin" 

C and the boundary of each panel is drawn with color "klredg". 

C NOTE: The two routines use different hidden line removal algorithm. 
C Which you use is a matter of taste. However the equivalent 

C of Q3DFIL is not available on the pen plotter since it is 

e based on an overpainting scheme which isn't practical on a plotter. 
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The angles phi and theta specify the view direction as shown in the 
appropriate figure in your GRAFMATIC/PLOTMATIC manual. 
The parameter icross is one if you wish cross-grid lines in both the 
x and y directions or zero for grid lines only parallel to 
the x axis (as may be approrpiate in creating a spectral plot). 


16 WRITE(*,*) ‘Enter 0 for 3D-stick or 1 for 3D-fill option’ 
READ(*,*) IFIL 
IFIL=0 
WRITE(*,*) 
WRITE(*,*) ' Default view angles are: phiz-45.,theta-70.' 
IF (IFIL.EQ.1) THEN 


(3009 S C)O). (6Y CY OI 


C KLRIN is the color inside the panel 
KLRIN=2 


C KLREDG is the panel border color 


E 
KLREDG=1 
ELSE 
C ICROSS is 1 for cross-grids in both directions 
C 
ICROSS=1 
ENDIF 


C Continuation point to loop through desired angles of observation 
C 
17 CONTINUE 
CALL QSMODE(MODEP) 


C Modify THETA to reflect YFIN/ZFIN 
€ 
STHETA = ATAN( (YFIN/ZFIN)*TAN(THETA*DTR) )/DTR 


C Q3DROT changes the view perspective as described in the documentation 
C 
CALL Q3DROT(X,Y,Z,M,N,PHI,STHETA) 
C 
CALL Q3DWIND(XST,XFIN,YST, YFIN,ZST,ZFIN,XMIN,XMAX,YMIN,YMAX) 
IOPT=0 
IF ((XMAX-XMIN).NE.0.0) THEN 
YOVERX = (YMAX-YMIN) / (XMAX-XMIN) 
ELSE 
YOVERX = 999 
ENDIF 
XORG=0.0 
YORG=0.0 
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CALL QPLOT(100,550,30,270,XMIN,XMAX, Y MIN, YMAX, 


1 XORG,YORG,IOPT, YOVERX,ASPECT) 


C The next 3 routines create the x, y, and z axes 
C if Q3DFIL is used, then draw the axes AFTER the 3D plot is drawn 
C so that the 3D plot does not overwrite the x/y axes. 


C 


C 


IF (IFIL.EQ.1) GO TO 66 


CALL Q3DXAX(XST,XFIN,XSTEP,9,1,2, YBEG, YEND,ZFLOOR, 1.0) 
CALL Q3DYAX(YST,YFIN, YSTEP,9,1,2,X BEG,XEND,ZFLOOR, 1.0) 
CALL Q3DZAX(ZST,ZFIN,ZSTEP,9,-1,2,XBEG, YBEG, 1.0) 


66 CONTINUE 


C Finally draw the 3D plot, using the appropriate option. 


t 


IF(IFIL.EQ.1) CALL Q3DFIL(X,Y,M,N,KLRIN,KLREDG) 
IF(IFIL.NE.1) CALL Q3DSTk(X,Y,M,N,IWORK1,IWORK2,NPTS,ICROSS) 


C Draw the axes if using Q3DFIL routine 


E 


C 
C 


IF(IFIL.NE.1) GO TO 67 

CALL Q3DXAX(XST,XFIN,XSTEP9,1,2, YBEG,YEND,ZFLOOR, 1.0) 
CALL Q3DYAX(YST, YFIN, YSTEP.9,1,2,XBEG,XEND,ZFLOOR, 1.0) 
CALL Q3DZAX(ZST,ZFIN,ZSTEP,9,-1,2,XBEG, YBEG, 1.0) 


67 CONTINUE 


2 GG GO GGG GO OG ak k ak ak ak ak a ak ak ak 3k ak k ak ak k ak Ok ak k ak k k 3k k k k k k k 


WRITE TEXT ON SCREEN 

NCHAR = 64 

KOLOR = 10 

ICOL = 10 

IROW = 24 

CALL QPTXT(NCHAR,TITLE 1,KOLOR,ICOL,IROW) 
SOO ke ok jk ok oko aooaa oaao ooo eee eek kk eee doi icici ak 
WRITE TEXT ON SCREEN 

NCHAR = 64 

KOLOR = 10 

ICOL = 10 

IROW = 23 

CALL QPTXT(NCHAR,HDR1,KOLOR,ICOL,IROW) 
SG GGG GOI Iai Iai ii ial iiaigigiai ai aigioiaiak ak ak ak iakak ak kak kak ak ak kakaiak kak kkk k k 
WRITE TEXT ON SCREEN 

NCHAR = 64 

KOLOR = 10 

ICOL = 10 

IROW = 22 

CALL QPTXT(NCHAR,HDR2,KOLOR,ICOL,IROW) 


ak ak ak ak kc ae kc ak ke he hee he Ie he hee he he hee We he he Fk okc oe ke ok ok okcokc ok oc oc ok He ke oc ke oe kc ok ke a ole a ae e ke oc oe oie kc okcok oc oe oie ok oc oe ol al a 


WRITE TEXT ON SCREEN 


1959 


00O 


(I CIC) 


NCHAR = 64 

KOLOR = 10 

ICOL = 10 

IROW =21 

CALL QPTXT(NCHAR,GRFLAB,KOLOR,ICOL,IROW) 
sk jc sj je je je jeje dc kc eek je jk oe oc kc jc ac je oe jc je oc kc c kc ccc eek ake ake ake a ak ake 3 ake ak ake ake a ake ak ake ak ake ake ak kc kc ec ke eco ok 
WRITE TEXT ON SCREEN 

NCHAR = 64 

KOLOR = 10 

ICOL = 15 

IROW = 20 

CALL QPTXT(NCHAR, TITLE2,KOLOR,ICOL,IROW) 


ak ak ak ak ak ak ale ae ade al al ae ade ade de ae ade ae ade ade ale ae ade ade al ad al ae e ad Ak al al e e ae e ad e ae ade ae e ae e al a a a ae e al a a ae a a a a a k k k k k k 


MAGNITUDES INTO A STRING 


CALL NUMSTR(ZMIN,MIN) 
CALL NUMSTR(ZMAX,MAX) 
CALL NUMSTR(ZAVG,AVG) 
TITLE3='min. Z = //MIN//; max. Z 2 //MAX 
^ /f; average Z 2 //AVG 

WRITE TEXT ON SCREEN 


NCHAR = 64 

KOLOR = 10 

ICOL = 17 

IROW = 3 

CALL QPTXT(NCHAR,TITLE3,KOLOR,ICOL,IROW) 


ak ak ade ade ade ale ale ade ale ale ade ale ade ade ale ade ale ade ale ade ade ade ale ale ade ade ale ale ale ale ade ade ale ale ale ale ale ade ale ade ale ale ale ade ade ale ale ade ale ak ade ade ale ade ale ie ade ade ae ade ade ale ade ale ale x 


ANGLES INTO A STRING 


CALL NUMSTR(ALPHA,ASTRING) 

CALL NUMSTR(GAMMA,GSTRING) 

CALL NUMSTR(THETA,TSTRING) 

CALL NUMSTR(PHI,PSTRING) 

TITLE4="Inc. A = '//ASTRING//;Plane G = '//GSTRING 

+ /f; View: phi=//PSTRING// theta=//TSTRING 
WRITE TEXT ON SCREEN 


NCHAR = 64 

KOLOR = 10 

ICOL = 10 

IROW = 2 

CALL QPTXT(NCHAR,TITLE4,KOLOR,ICOL,IROW) 
al ad ade ade ade ale ad ade al al ade ad ad ad al al al ad ad al al ad a a ad al a al al a e ad ad ad ad al al a a al a a ad ad al al al a al a a a a ad ad al al al al ate ae de al e le al 
Hardcopy Query 

NCHAR = 30 

KOLOR = 12 

ICOL = 1 

IROW = 1 

HCOPY='Hardcopy ---> Enter P or p' 

CALL QPTXT(NCHAR,HCOPY KOLOR,ICOL,IROW) 

CALL QCMOV (32,1) 
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CALL QINKEY(IEXTEND,IKEY) 

NCHAR = 40 

KOLOR = 1 

ICOL = 1 

IROW = 1 

HCOPY=' 

CALL QPTXT(NCHAR,HCOPY,KOLOR,ICOL,IROW) 
C If 'P'or 'p' then prtsc 


IF( IKEY.EQ.80 .OR. IKEY.EQ.112 ) THEN 
CALL QPSCRN 
ENDIF 
C ak ak afe afe oic ok ok ofc ok kc ok aic oc ok oic ofc oc ok oic aic ok ok aic aic ok aic ok aic oc ofc oic ok oc ic ok fc oc ok afe afe aic aic ok ok ok ofc oic aic ok ofc ok kc ok ofc ok ok ake afe ok ofc okc ok okc ok oco 

CALL Q3DINV(X,Y,Z,M,N) 
CALL QSMODE(2) 
WRITE(*,*)'"ENTER NEW ANGLES PHI & THETA (DEG) OR (9,9) WHEN DONE 
READ(5,*)PHI,THETA 
IF ( (PHI.NE.9.).AND.(THETA.NE.9.) ) GOTO 17 

999 CONTINUE 

100 FORMAT(A1) 

520 FORMAT(3X,I5) 

578 FORMAT(3X,E12.6) 


101 FORMAT(A) 
RETURN 
END 


ak ak aak ak 3k ak i i Gi al a a ak ak ak ak 


E e ode al ad aj al aj a ad al a oc oic oc ode ad a a aj a aj a a a a a al a ad ad a al a a ad a a oic oic oc ok a a a a ad ad al a a a ad ae ae fe a ae ad al a a 


C 
SUBROUTINE XYSET(X,Y,M,N) 
c 
REAL X(MN),Y(M,N),I1,J1,N2,M2 


E 
C define x,y values over the m by n grid 
C 


M2-FLOAT(M-1)/2.0 
N2-FLOAT(N-1)/2.0 
DO 111 I=1,M 
DO 121 J=1,N 
11=FLOAT(D-1.0 
J1=FLOAT()-1.0 
X(1J)= (11-M2)/M2 
Y (1,)= (J1-N2)/N2 
121 CONTINUE 
111 CONTINUE 
RETURN 
END 
ILLO E E E al E e A e e e E e ae ak ake aka NE 
C 
C SUBROUTINE NUMSTR 
E 
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C ACCEPTS NUMBER AND RETURNS CHAR*5 WITH SIGN AND THREE DIGIT STRING 


SUBROUTINE NUMSTR(ANUM,STR) 
REAL ANUM,VNUM 

INTEGER NUM,HUN,TEN,ONE,ORDER 
CHARACTER*5 STR 


ORDER = 0 
C Get sign information 


IF (ANUM.GE.0.0) THEN 
STR(1:1)='+' 

ELSE 
STR(1:1)='-' 

ENDIF 


C Strip sign values 
VNUM=ABS(ANUM) 


IF (VNUM .LT. 1.0) THEN 
ORDER=0 
VNUM=VNUM* 1000.0 

ELSEIF (VNUM .LT .10.0) THEN 
ORDER-1 

VNUM=VNUM* 100.0 

ELSEIF (VNUM .LT. 100.0) THEN 
ORDER=2 
VNUM=VNUM*10.0 

ELSE 

ORDER=3 

ENDIF 


C Use integer part of shifted VNUM 


NUM=VNUM 

IF (NUM.GE.1000) THEN NUM=999 
HUN=INT(NUM/100.) 
TEN=(NUM-100* HUN)/10 
ONE=(NUM-100*HUN-10* TEN) 


C PLACE DECIMAL 

C 
IF (ORDER .EQ. 0) THEN 
STR(2:2)=".' 
STR(3:3)=CHAR(48+HUN) 
STR(4:4)=CHAR(48+TEN) 
STR(5:5)=CHAR(48+ONE) 
ELSEIF (ORDER .EQ. 1) THEN 
STR(2:2)=CHAR(48+HUN) 
SERO 3)=. 
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STR(4:4)=CHAR(48+TEN) 
STR(S:5)=CHAR(48+0NE) 
ELSEIF (ORDER .EQ. 2) THEN 
STR(2:2)=CHAR(48+HUN) 
STR(3:3)=CHAR(48+TEN) 
STR(4:4)=".' 
STR(S:5)=CHAR(48+0NE) 
ELSE 
STR(2:2)=CHAR(48+HUN) 
STR(3:3)=CHAR(48+TEN) 
STR(4:4)=CHAR(48+0NE) 
STR(5:5)=".' 
ENDIF 


RETURN 
END 
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Appendix D 


SOFTWARE SOURCES 
Microsoft FORTRAN 
16011 NE 36th Way 
Box 97017 


Redmond, WA 98073-9717 


West Coast Consultants CURVE DIGITIZER 
4202 Genesee Avenue, Suite 309 
San Diego, CA 92117 
(619) 565-1266 


Microcompatibles, Inc. GRAFMATIC 
301 Prelude Drive 
Silver Spring, MD 20901 
(301)593-5151 


Prof. M.A. Morgan 
Code 62Mw 
Naval Postgraduate School 
Monterey, CA 93943 
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